13 int libsvm_version = LIBSVM_VERSION;
15 typedef signed char schar;
17 template <
class T>
static inline T min(T x,T y) {
return (x<y)?x:y; }
20 template <
class T>
static inline T max(T x,T y) {
return (x>y)?x:y; }
22 template <
class T>
static inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
23 template <
class S,
class T>
static inline void clone(T*& dst, S* src,
int n)
26 memcpy((
void *)dst,(
void *)src,
sizeof(T)*n);
28 static inline double powi(
double base,
int times)
30 double tmp = base, ret = 1.0;
32 for(
int t=times; t>0; t/=2)
41 #define Malloc(type,n) (type *)malloc((n)*sizeof(type)) 43 static void print_string_stdout(
const char *s)
48 static void (*svm_print_string) (
const char *) = &print_string_stdout;
50 static void info(
const char *fmt,...)
57 (*svm_print_string)(buf);
60 static void info(
const char *fmt,...) {}
72 Cache(
int l,
long int size);
78 int get_data(
const int index, Qfloat **data,
int len);
79 void swap_index(
int i,
int j);
92 void lru_delete(head_t *h);
93 void lru_insert(head_t *h);
96 Cache::Cache(
int l_,
long int size_):l(l_),size(size_)
98 head = (head_t *)calloc(l,
sizeof(head_t));
99 size /=
sizeof(Qfloat);
100 size -= l *
sizeof(head_t) /
sizeof(Qfloat);
101 size = max(size, 2 * (
long int) l);
102 lru_head.next = lru_head.prev = &lru_head;
107 for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
112 void Cache::lru_delete(head_t *h)
115 h->prev->next = h->next;
116 h->next->prev = h->prev;
119 void Cache::lru_insert(head_t *h)
123 h->prev = lru_head.prev;
128 int Cache::get_data(
const int index, Qfloat **data,
int len)
130 head_t *h = &head[index];
131 if(h->len) lru_delete(h);
132 int more = len - h->len;
139 head_t *old = lru_head.next;
148 h->data = (Qfloat *)realloc(h->data,
sizeof(Qfloat)*len);
158 void Cache::swap_index(
int i,
int j)
162 if(head[i].len) lru_delete(&head[i]);
163 if(head[j].len) lru_delete(&head[j]);
164 swap(head[i].data,head[j].data);
165 swap(head[i].len,head[j].len);
166 if(head[i].len) lru_insert(&head[i]);
167 if(head[j].len) lru_insert(&head[j]);
170 for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
175 swap(h->data[i],h->data[j]);
198 virtual Qfloat *get_Q(
int column,
int len)
const = 0;
199 virtual double *get_QD()
const = 0;
200 virtual void swap_index(
int i,
int j)
const = 0;
211 virtual Qfloat *get_Q(
int column,
int len)
const = 0;
212 virtual double *get_QD()
const = 0;
213 virtual void swap_index(
int i,
int j)
const 216 if(x_square) swap(x_square[i],x_square[j]);
220 double (
Kernel::*kernel_function)(
int i,
int j)
const;
227 const int kernel_type;
233 double kernel_linear(
int i,
int j)
const 235 return dot(x[i],x[j]);
237 double kernel_poly(
int i,
int j)
const 239 return powi(gamma*dot(x[i],x[j])+coef0,degree);
241 double kernel_rbf(
int i,
int j)
const 243 return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
245 double kernel_sigmoid(
int i,
int j)
const 247 return tanh(gamma*dot(x[i],x[j])+coef0);
249 double kernel_precomputed(
int i,
int j)
const 251 return x[i][(int)(x[j][0].value)].value;
256 :kernel_type(param.kernel_type), degree(param.degree),
257 gamma(param.gamma), coef0(param.coef0)
262 kernel_function = &Kernel::kernel_linear;
265 kernel_function = &Kernel::kernel_poly;
268 kernel_function = &Kernel::kernel_rbf;
271 kernel_function = &Kernel::kernel_sigmoid;
274 kernel_function = &Kernel::kernel_precomputed;
280 if(kernel_type == RBF)
282 x_square =
new double[l];
284 x_square[i] = dot(x[i],x[i]);
299 while(px->index != -1 && py->index != -1)
301 if(px->index == py->index)
303 sum += px->value * py->value;
309 if(px->index > py->index)
321 switch(param.kernel_type)
326 return powi(param.gamma*dot(x,y)+param.coef0,param.degree);
330 while(x->index != -1 && y->index !=-1)
332 if(x->index == y->index)
334 double d = x->value - y->value;
341 if(x->index > y->index)
343 sum += y->value * y->value;
348 sum += x->value * x->value;
354 while(x->index != -1)
356 sum += x->value * x->value;
360 while(y->index != -1)
362 sum += y->value * y->value;
366 return exp(-param.gamma*sum);
369 return tanh(param.gamma*dot(x,y)+param.coef0);
371 return x[(int)(y->value)].value;
403 double upper_bound_p;
404 double upper_bound_n;
408 void Solve(
int l,
const QMatrix& Q,
const double *p_,
const schar *y_,
409 double *alpha_,
double Cp,
double Cn,
double eps,
415 enum { LOWER_BOUND, UPPER_BOUND, FREE };
430 return (y[i] > 0)? Cp : Cn;
432 void update_alpha_status(
int i)
434 if(alpha[i] >= get_C(i))
435 alpha_status[i] = UPPER_BOUND;
436 else if(alpha[i] <= 0)
437 alpha_status[i] = LOWER_BOUND;
438 else alpha_status[i] = FREE;
440 bool is_upper_bound(
int i) {
return alpha_status[i] == UPPER_BOUND; }
441 bool is_lower_bound(
int i) {
return alpha_status[i] == LOWER_BOUND; }
442 bool is_free(
int i) {
return alpha_status[i] == FREE; }
443 void swap_index(
int i,
int j);
444 void reconstruct_gradient();
445 virtual int select_working_set(
int &i,
int &j);
446 virtual double calculate_rho();
447 virtual void do_shrinking();
449 bool be_shrunk(
int i,
double Gmax1,
double Gmax2);
452 void Solver::swap_index(
int i,
int j)
457 swap(alpha_status[i],alpha_status[j]);
458 swap(alpha[i],alpha[j]);
460 swap(active_set[i],active_set[j]);
461 swap(G_bar[i],G_bar[j]);
464 void Solver::reconstruct_gradient()
468 if(active_size == l)
return;
473 for(j=active_size;j<l;j++)
474 G[j] = G_bar[j] + p[j];
476 for(j=0;j<active_size;j++)
480 if(2*nr_free < active_size)
481 info(
"\nWARNING: using -h 0 may be faster\n");
483 if (nr_free*l > 2*active_size*(l-active_size))
485 for(i=active_size;i<l;i++)
487 const Qfloat *Q_i = Q->get_Q(i,active_size);
488 for(j=0;j<active_size;j++)
490 G[i] += alpha[j] * Q_i[j];
495 for(i=0;i<active_size;i++)
498 const Qfloat *Q_i = Q->get_Q(i,l);
499 double alpha_i = alpha[i];
500 for(j=active_size;j<l;j++)
501 G[j] += alpha_i * Q_i[j];
506 void Solver::Solve(
int l,
const QMatrix& Q,
const double *p_,
const schar *y_,
507 double *alpha_,
double Cp,
double Cn,
double eps,
516 clone(alpha,alpha_,l);
524 alpha_status =
new char[l];
526 update_alpha_status(i);
531 active_set =
new int[l];
540 G_bar =
new double[l];
548 if(!is_lower_bound(i))
550 const Qfloat *Q_i = Q.get_Q(i,l);
551 double alpha_i = alpha[i];
554 G[j] += alpha_i*Q_i[j];
555 if(is_upper_bound(i))
557 G_bar[j] += get_C(i) * Q_i[j];
564 int max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);
565 int counter = min(l,1000)+1;
567 while(iter < max_iter)
573 counter = min(l,1000);
574 if(shrinking) do_shrinking();
580 if(select_working_set(i,j)!=0)
583 reconstruct_gradient();
588 if(select_working_set(i,j)!=0)
598 const Qfloat *Q_i = Q.get_Q(i,active_size);
599 const Qfloat *Q_j = Q.get_Q(j,active_size);
601 double C_i = get_C(i);
602 double C_j = get_C(j);
604 double old_alpha_i = alpha[i];
605 double old_alpha_j = alpha[j];
609 double quad_coef = QD[i]+QD[j]+2*Q_i[j];
612 double delta = (-G[i]-G[j])/quad_coef;
613 double diff = alpha[i] - alpha[j];
638 alpha[j] = C_i - diff;
646 alpha[i] = C_j + diff;
652 double quad_coef = QD[i]+QD[j]-2*Q_i[j];
655 double delta = (G[i]-G[j])/quad_coef;
656 double sum = alpha[i] + alpha[j];
665 alpha[j] = sum - C_i;
681 alpha[i] = sum - C_j;
696 double delta_alpha_i = alpha[i] - old_alpha_i;
697 double delta_alpha_j = alpha[j] - old_alpha_j;
699 for(
int k=0;k<active_size;k++)
701 G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
707 bool ui = is_upper_bound(i);
708 bool uj = is_upper_bound(j);
709 update_alpha_status(i);
710 update_alpha_status(j);
712 if(ui != is_upper_bound(i))
717 G_bar[k] -= C_i * Q_i[k];
720 G_bar[k] += C_i * Q_i[k];
723 if(uj != is_upper_bound(j))
728 G_bar[k] -= C_j * Q_j[k];
731 G_bar[k] += C_j * Q_j[k];
741 reconstruct_gradient();
746 info(
"\nWARNING: reaching max number of iterations");
751 si->rho = calculate_rho();
758 v += alpha[i] * (G[i] + p[i]);
766 alpha_[active_set[i]] = alpha[i];
777 si->upper_bound_p = Cp;
778 si->upper_bound_n = Cn;
781 info(
"\noptimization finished, #iter = %d\n",iter);
786 delete[] alpha_status;
793 int Solver::select_working_set(
int &out_i,
int &out_j)
805 double obj_diff_min = INF;
807 for(
int t=0;t<active_size;t++)
810 if(!is_upper_bound(t))
819 if(!is_lower_bound(t))
828 const Qfloat *Q_i = NULL;
830 Q_i = Q->get_Q(i,active_size);
832 for(
int j=0;j<active_size;j++)
836 if (!is_lower_bound(j))
838 double grad_diff=Gmax+G[j];
844 double quad_coef = QD[i]+QD[j]-2.0*y[i]*Q_i[j];
846 obj_diff = -(grad_diff*grad_diff)/quad_coef;
848 obj_diff = -(grad_diff*grad_diff)/TAU;
850 if (obj_diff <= obj_diff_min)
853 obj_diff_min = obj_diff;
860 if (!is_upper_bound(j))
862 double grad_diff= Gmax-G[j];
868 double quad_coef = QD[i]+QD[j]+2.0*y[i]*Q_i[j];
870 obj_diff = -(grad_diff*grad_diff)/quad_coef;
872 obj_diff = -(grad_diff*grad_diff)/TAU;
874 if (obj_diff <= obj_diff_min)
877 obj_diff_min = obj_diff;
892 bool Solver::be_shrunk(
int i,
double Gmax1,
double Gmax2)
894 if(is_upper_bound(i))
897 return(-G[i] > Gmax1);
899 return(-G[i] > Gmax2);
901 else if(is_lower_bound(i))
904 return(G[i] > Gmax2);
906 return(G[i] > Gmax1);
912 void Solver::do_shrinking()
919 for(i=0;i<active_size;i++)
923 if(!is_upper_bound(i))
928 if(!is_lower_bound(i))
936 if(!is_upper_bound(i))
941 if(!is_lower_bound(i))
949 if(unshrink ==
false && Gmax1 + Gmax2 <= eps*10)
952 reconstruct_gradient();
957 for(i=0;i<active_size;i++)
958 if (be_shrunk(i, Gmax1, Gmax2))
961 while (active_size > i)
963 if (!be_shrunk(active_size, Gmax1, Gmax2))
965 swap_index(i,active_size);
973 double Solver::calculate_rho()
977 double ub = INF, lb = -INF, sum_free = 0;
978 for(
int i=0;i<active_size;i++)
980 double yG = y[i]*G[i];
982 if(is_upper_bound(i))
989 else if(is_lower_bound(i))
1004 r = sum_free/nr_free;
1020 void Solve(
int l,
const QMatrix& Q,
const double *p,
const schar *y,
1021 double *alpha,
double Cp,
double Cn,
double eps,
1025 Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking,verbose);
1029 int select_working_set(
int &i,
int &j);
1030 double calculate_rho();
1031 bool be_shrunk(
int i,
double Gmax1,
double Gmax2,
double Gmax3,
double Gmax4);
1032 void do_shrinking();
1036 int Solver_NU::select_working_set(
int &out_i,
int &out_j)
1044 double Gmaxp = -INF;
1045 double Gmaxp2 = -INF;
1048 double Gmaxn = -INF;
1049 double Gmaxn2 = -INF;
1053 double obj_diff_min = INF;
1055 for(
int t=0;t<active_size;t++)
1058 if(!is_upper_bound(t))
1067 if(!is_lower_bound(t))
1077 const Qfloat *Q_ip = NULL;
1078 const Qfloat *Q_in = NULL;
1080 Q_ip = Q->get_Q(ip,active_size);
1082 Q_in = Q->get_Q(in,active_size);
1084 for(
int j=0;j<active_size;j++)
1088 if (!is_lower_bound(j))
1090 double grad_diff=Gmaxp+G[j];
1096 double quad_coef = QD[ip]+QD[j]-2*Q_ip[j];
1098 obj_diff = -(grad_diff*grad_diff)/quad_coef;
1100 obj_diff = -(grad_diff*grad_diff)/TAU;
1102 if (obj_diff <= obj_diff_min)
1105 obj_diff_min = obj_diff;
1112 if (!is_upper_bound(j))
1114 double grad_diff=Gmaxn-G[j];
1115 if (-G[j] >= Gmaxn2)
1120 double quad_coef = QD[in]+QD[j]-2*Q_in[j];
1122 obj_diff = -(grad_diff*grad_diff)/quad_coef;
1124 obj_diff = -(grad_diff*grad_diff)/TAU;
1126 if (obj_diff <= obj_diff_min)
1129 obj_diff_min = obj_diff;
1136 if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
1139 if (y[Gmin_idx] == +1)
1148 bool Solver_NU::be_shrunk(
int i,
double Gmax1,
double Gmax2,
double Gmax3,
double Gmax4)
1150 if(is_upper_bound(i))
1153 return(-G[i] > Gmax1);
1155 return(-G[i] > Gmax4);
1157 else if(is_lower_bound(i))
1160 return(G[i] > Gmax2);
1162 return(G[i] > Gmax3);
1168 void Solver_NU::do_shrinking()
1170 double Gmax1 = -INF;
1171 double Gmax2 = -INF;
1172 double Gmax3 = -INF;
1173 double Gmax4 = -INF;
1177 for(i=0;i<active_size;i++)
1179 if(!is_upper_bound(i))
1183 if(-G[i] > Gmax1) Gmax1 = -G[i];
1185 else if(-G[i] > Gmax4) Gmax4 = -G[i];
1187 if(!is_lower_bound(i))
1191 if(G[i] > Gmax2) Gmax2 = G[i];
1193 else if(G[i] > Gmax3) Gmax3 = G[i];
1197 if(unshrink ==
false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10)
1200 reconstruct_gradient();
1204 for(i=0;i<active_size;i++)
1205 if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
1208 while (active_size > i)
1210 if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
1212 swap_index(i,active_size);
1220 double Solver_NU::calculate_rho()
1222 int nr_free1 = 0,nr_free2 = 0;
1223 double ub1 = INF, ub2 = INF;
1224 double lb1 = -INF, lb2 = -INF;
1225 double sum_free1 = 0, sum_free2 = 0;
1227 for(
int i=0;i<active_size;i++)
1231 if(is_upper_bound(i))
1232 lb1 = max(lb1,G[i]);
1233 else if(is_lower_bound(i))
1234 ub1 = min(ub1,G[i]);
1243 if(is_upper_bound(i))
1244 lb2 = max(lb2,G[i]);
1245 else if(is_lower_bound(i))
1246 ub2 = min(ub2,G[i]);
1257 r1 = sum_free1/nr_free1;
1262 r2 = sum_free2/nr_free2;
1277 :
Kernel(prob.l, prob.x, param)
1280 cache =
new Cache(prob.l,(
long int)(param.cache_size*(1<<20)));
1281 QD =
new double[prob.l];
1282 for(
int i=0;i<prob.l;i++)
1283 QD[i] = (this->*kernel_function)(i,i);
1286 Qfloat *get_Q(
int i,
int len)
const 1290 if((start = cache->get_data(i,&data,len)) < len)
1292 for(j=start;j<len;j++)
1293 data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
1298 double *get_QD()
const 1303 void swap_index(
int i,
int j)
const 1305 cache->swap_index(i,j);
1306 Kernel::swap_index(i,j);
1327 :
Kernel(prob.l, prob.x, param)
1329 cache =
new Cache(prob.l,(
long int)(param.cache_size*(1<<20)));
1330 QD =
new double[prob.l];
1331 for(
int i=0;i<prob.l;i++)
1332 QD[i] = (this->*kernel_function)(i,i);
1335 Qfloat *get_Q(
int i,
int len)
const 1339 if((start = cache->get_data(i,&data,len)) < len)
1341 for(j=start;j<len;j++)
1342 data[j] = (Qfloat)(this->*kernel_function)(i,j);
1347 double *get_QD()
const 1352 void swap_index(
int i,
int j)
const 1354 cache->swap_index(i,j);
1355 Kernel::swap_index(i,j);
1373 :
Kernel(prob.l, prob.x, param)
1376 cache =
new Cache(l,(
long int)(param.cache_size*(1<<20)));
1377 QD =
new double[2*l];
1378 sign =
new schar[2*l];
1379 index =
new int[2*l];
1380 for(
int k=0;k<l;k++)
1386 QD[k] = (this->*kernel_function)(k,k);
1389 buffer[0] =
new Qfloat[2*l];
1390 buffer[1] =
new Qfloat[2*l];
1394 void swap_index(
int i,
int j)
const 1396 swap(sign[i],sign[j]);
1397 swap(index[i],index[j]);
1401 Qfloat *get_Q(
int i,
int len)
const 1404 int j, real_i = index[i];
1405 if(cache->get_data(real_i,&data,l) < l)
1408 data[j] = (Qfloat)(this->*kernel_function)(real_i,j);
1412 Qfloat *buf = buffer[next_buffer];
1413 next_buffer = 1 - next_buffer;
1416 buf[j] = (Qfloat) si * (Qfloat) sign[j] * data[index[j]];
1420 double *get_QD()
const 1439 mutable int next_buffer;
1447 static void solve_c_svc(
1452 double *minus_ones =
new double[l];
1453 schar *y =
new schar[l];
1461 if(prob->y[i] > 0) y[i] = +1;
else y[i] = -1;
1465 s.Solve(l,
SVC_Q(*prob,*param,y), minus_ones, y,
1466 alpha, Cp, Cn, param->eps, si, param->shrinking, param->verbose);
1470 sum_alpha += alpha[i];
1474 info(
"nu = %f\n", sum_alpha/(Cp*prob->l));
1479 delete[] minus_ones;
1483 static void solve_nu_svc(
1489 double nu = param->nu;
1491 schar *y =
new schar[l];
1499 double sum_pos = nu*l/2;
1500 double sum_neg = nu*l/2;
1505 alpha[i] = min(1.0,sum_pos);
1506 sum_pos -= alpha[i];
1510 alpha[i] = min(1.0,sum_neg);
1511 sum_neg -= alpha[i];
1514 double *zeros =
new double[l];
1520 s.Solve(l,
SVC_Q(*prob,*param,y), zeros, y,
1521 alpha, 1.0, 1.0, param->eps, si, param->shrinking, param->verbose);
1525 info(
"C = %f\n",1/r);
1532 si->upper_bound_p = 1/r;
1533 si->upper_bound_n = 1/r;
1539 static void solve_one_class(
1544 double *zeros =
new double[l];
1545 schar *ones =
new schar[l];
1548 int n = (int)(param->nu*prob->l);
1553 alpha[n] = param->nu * prob->l - n;
1564 s.Solve(l,
ONE_CLASS_Q(*prob,*param), zeros, ones,
1565 alpha, 1.0, 1.0, param->eps, si, param->shrinking, param->verbose );
1571 static void solve_epsilon_svr(
1576 double *alpha2 =
new double[2*l];
1577 double *linear_term =
new double[2*l];
1578 schar *y =
new schar[2*l];
1584 linear_term[i] = param->p - prob->y[i];
1588 linear_term[i+l] = param->p + prob->y[i];
1593 s.Solve(2*l,
SVR_Q(*prob,*param), linear_term, y,
1594 alpha2, param->C, param->C, param->eps, si, param->shrinking, param->verbose);
1596 double sum_alpha = 0;
1599 alpha[i] = alpha2[i] - alpha2[i+l];
1600 sum_alpha += fabs(alpha[i]);
1603 info(
"nu = %f\n",sum_alpha/(param->C*l));
1606 delete[] linear_term;
1610 static void solve_nu_svr(
1615 double C = param->C;
1616 double *alpha2 =
new double[2*l];
1617 double *linear_term =
new double[2*l];
1618 schar *y =
new schar[2*l];
1621 double sum = C * param->nu * l / 2;
1624 alpha2[i] = alpha2[i+l] = min(sum,C);
1627 linear_term[i] = - prob->y[i];
1630 linear_term[i+l] = prob->y[i];
1635 s.Solve(2*l,
SVR_Q(*prob,*param), linear_term, y,
1636 alpha2, C, C, param->eps, si, param->shrinking, param->verbose);
1639 info(
"epsilon = %f\n",-si->r);
1642 alpha[i] = alpha2[i] - alpha2[i+l];
1645 delete[] linear_term;
1660 double Cp,
double Cn)
1662 double *alpha = Malloc(
double,prob->l);
1664 switch(param->svm_type)
1667 solve_c_svc(prob,param,alpha,&si,Cp,Cn);
1670 solve_nu_svc(prob,param,alpha,&si);
1673 solve_one_class(prob,param,alpha,&si);
1676 solve_epsilon_svr(prob,param,alpha,&si);
1679 solve_nu_svr(prob,param,alpha,&si);
1684 info(
"obj = %f, rho = %f\n",si.obj,si.rho);
1690 for(
int i=0;i<prob->l;i++)
1692 if(fabs(alpha[i]) > 0)
1697 if(fabs(alpha[i]) >= si.upper_bound_p)
1702 if(fabs(alpha[i]) >= si.upper_bound_n)
1709 info(
"nSV = %d, nBSV = %d\n",nSV,nBSV);
1718 static void sigmoid_train(
1719 int l,
const double *dec_values,
const double *labels,
1720 double& A,
double& B)
1722 double prior1=0, prior0 = 0;
1726 if (labels[i] > 0) prior1+=1;
1730 double min_step=1e-10;
1733 double hiTarget=(prior1+1.0)/(prior1+2.0);
1734 double loTarget=1/(prior0+2.0);
1735 double *t=Malloc(
double,l);
1736 double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
1737 double newA,newB,newf,d1,d2;
1741 A=0.0; B=log((prior0+1.0)/(prior1+1.0));
1746 if (labels[i]>0) t[i]=hiTarget;
1748 fApB = dec_values[i]*A+B;
1750 fval += t[i]*fApB + log(1+exp(-fApB));
1752 fval += (t[i] - 1)*fApB +log(1+exp(fApB));
1754 for (iter=0;iter<max_iter;iter++)
1759 h21=0.0;g1=0.0;g2=0.0;
1762 fApB = dec_values[i]*A+B;
1765 p=exp(-fApB)/(1.0+exp(-fApB));
1766 q=1.0/(1.0+exp(-fApB));
1770 p=1.0/(1.0+exp(fApB));
1771 q=exp(fApB)/(1.0+exp(fApB));
1774 h11+=dec_values[i]*dec_values[i]*d2;
1776 h21+=dec_values[i]*d2;
1778 g1+=dec_values[i]*d1;
1783 if (fabs(g1)<eps && fabs(g2)<eps)
1787 det=h11*h22-h21*h21;
1788 dA=-(h22*g1 - h21 * g2) / det;
1789 dB=-(-h21*g1+ h11 * g2) / det;
1794 while (stepsize >= min_step)
1796 newA = A + stepsize * dA;
1797 newB = B + stepsize * dB;
1803 fApB = dec_values[i]*newA+newB;
1805 newf += t[i]*fApB + log(1+exp(-fApB));
1807 newf += (t[i] - 1)*fApB +log(1+exp(fApB));
1810 if (newf<fval+0.0001*stepsize*gd)
1812 A=newA;B=newB;fval=newf;
1816 stepsize = stepsize / 2.0;
1819 if (stepsize < min_step)
1821 info(
"Line search fails in two-class probability estimates\n");
1827 info(
"Reaching maximal iterations in two-class probability estimates\n");
1831 static double sigmoid_predict(
double decision_value,
double A,
double B)
1833 double fApB = decision_value*A+B;
1836 return exp(-fApB)/(1.0+exp(-fApB));
1838 return 1.0/(1+exp(fApB)) ;
1842 static void multiclass_probability(
int k,
double **r,
double *p)
1845 int iter = 0, max_iter=max(100,k);
1846 double **Q=Malloc(
double *,k);
1847 double *Qp=Malloc(
double,k);
1848 double pQp, eps=0.005/k;
1853 Q[t]=Malloc(
double,k);
1857 Q[t][t]+=r[j][t]*r[j][t];
1862 Q[t][t]+=r[j][t]*r[j][t];
1863 Q[t][j]=-r[j][t]*r[t][j];
1866 for (iter=0;iter<max_iter;iter++)
1874 Qp[t]+=Q[t][j]*p[j];
1880 double error=fabs(Qp[t]-pQp);
1881 if (error>max_error)
1884 if (max_error<eps)
break;
1888 double diff=(-Qp[t]+pQp)/Q[t][t];
1890 pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
1893 Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
1899 info(
"Exceeds max_iter in multiclass_prob\n");
1900 for(t=0;t<k;t++) free(Q[t]);
1906 static void svm_binary_svc_probability(
1908 double Cp,
double Cn,
double& probA,
double& probB)
1912 int *perm = Malloc(
int,prob->l);
1913 double *dec_values = Malloc(
double,prob->l);
1916 for(i=0;i<prob->l;i++) perm[i]=i;
1917 for(i=0;i<prob->l;i++)
1919 int j = i+rand()%(prob->l-i);
1920 swap(perm[i],perm[j]);
1922 for(i=0;i<nr_fold;i++)
1924 int begin = i*prob->l/nr_fold;
1925 int end = (i+1)*prob->l/nr_fold;
1929 subprob.l = prob->l-(end-begin);
1930 subprob.x = Malloc(
struct svm_node*,subprob.l);
1931 subprob.y = Malloc(
double,subprob.l);
1934 for(j=0;j<begin;j++)
1936 subprob.x[k] = prob->x[perm[j]];
1937 subprob.y[k] = prob->y[perm[j]];
1940 for(j=end;j<prob->l;j++)
1942 subprob.x[k] = prob->x[perm[j]];
1943 subprob.y[k] = prob->y[perm[j]];
1946 int p_count=0,n_count=0;
1953 if(p_count==0 && n_count==0)
1954 for(j=begin;j<end;j++)
1955 dec_values[perm[j]] = 0;
1956 else if(p_count > 0 && n_count == 0)
1957 for(j=begin;j<end;j++)
1958 dec_values[perm[j]] = 1;
1959 else if(p_count == 0 && n_count > 0)
1960 for(j=begin;j<end;j++)
1961 dec_values[perm[j]] = -1;
1965 subparam.probability=0;
1967 subparam.nr_weight=2;
1968 subparam.weight_label = Malloc(
int,2);
1969 subparam.weight = Malloc(
double,2);
1970 subparam.weight_label[0]=+1;
1971 subparam.weight_label[1]=-1;
1972 subparam.weight[0]=Cp;
1973 subparam.weight[1]=Cn;
1974 struct svm_model *submodel = svm_train(&subprob,&subparam);
1975 for(j=begin;j<end;j++)
1977 svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]]));
1979 dec_values[perm[j]] *= submodel->label[0];
1981 svm_free_and_destroy_model(&submodel);
1982 svm_destroy_param(&subparam);
1987 sigmoid_train(prob->l,dec_values,prob->y,probA,probB);
1993 static double svm_svr_probability(
1998 double *ymv = Malloc(
double,prob->l);
2002 newparam.probability = 0;
2003 svm_cross_validation(prob,&newparam,nr_fold,ymv);
2004 for(i=0;i<prob->l;i++)
2006 ymv[i]=prob->y[i]-ymv[i];
2007 mae += fabs(ymv[i]);
2010 double std=sqrt(2*mae*mae);
2013 for(i=0;i<prob->l;i++)
2014 if (fabs(ymv[i]) > 5*std)
2018 mae /= (prob->l-count);
2019 info(
"Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
2027 static void svm_group_classes(
const svm_problem *prob,
int *nr_class_ret,
int **label_ret,
int **start_ret,
int **count_ret,
int *perm)
2030 int max_nr_class = 16;
2032 int *label = Malloc(
int,max_nr_class);
2033 int *count = Malloc(
int,max_nr_class);
2034 int *data_label = Malloc(
int,l);
2039 int this_label = (int)prob->y[i];
2041 for(j=0;j<nr_class;j++)
2043 if(this_label == label[j])
2052 if(nr_class == max_nr_class)
2055 label = (
int *)realloc(label,max_nr_class*
sizeof(
int));
2056 count = (
int *)realloc(count,max_nr_class*
sizeof(
int));
2058 label[nr_class] = this_label;
2059 count[nr_class] = 1;
2064 int *start = Malloc(
int,nr_class);
2066 for(i=1;i<nr_class;i++)
2067 start[i] = start[i-1]+count[i-1];
2070 perm[start[data_label[i]]] = i;
2071 ++start[data_label[i]];
2074 for(i=1;i<nr_class;i++)
2075 start[i] = start[i-1]+count[i-1];
2077 *nr_class_ret = nr_class;
2090 model->param = *param;
2093 if(param->svm_type == ONE_CLASS ||
2094 param->svm_type == EPSILON_SVR ||
2095 param->svm_type == NU_SVR)
2098 model->nr_class = 2;
2099 model->label = NULL;
2101 model->probA = NULL; model->probB = NULL;
2102 model->sv_coef = Malloc(
double *,1);
2104 if(param->probability &&
2105 (param->svm_type == EPSILON_SVR ||
2106 param->svm_type == NU_SVR))
2108 model->probA = Malloc(
double,1);
2109 model->probA[0] = svm_svr_probability(prob,param);
2113 model->rho = Malloc(
double,1);
2114 model->rho[0] = f.rho;
2118 for(i=0;i<prob->l;i++)
2119 if(fabs(f.alpha[i]) > 0) ++nSV;
2121 model->SV = Malloc(
svm_node *,nSV);
2122 model->sv_coef[0] = Malloc(
double,nSV);
2124 for(i=0;i<prob->l;i++)
2125 if(fabs(f.alpha[i]) > 0)
2127 model->SV[j] = prob->x[i];
2128 model->sv_coef[0][j] = f.alpha[i];
2142 int *perm = Malloc(
int,l);
2145 svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
2147 info(
"WARNING: training data in only one class. See README for details.\n");
2152 x[i] = prob->x[perm[i]];
2156 double *weighted_C = Malloc(
double, nr_class);
2157 for(i=0;i<nr_class;i++)
2158 weighted_C[i] = param->C;
2159 for(i=0;i<param->nr_weight;i++)
2162 for(j=0;j<nr_class;j++)
2163 if(param->weight_label[i] == label[j])
2166 fprintf(stderr,
"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]);
2168 weighted_C[j] *= param->weight[i];
2173 bool *nonzero = Malloc(
bool,l);
2178 double *probA=NULL,*probB=NULL;
2179 if (param->probability)
2181 probA=Malloc(
double,nr_class*(nr_class-1)/2);
2182 probB=Malloc(
double,nr_class*(nr_class-1)/2);
2186 for(i=0;i<nr_class;i++)
2187 for(
int j=i+1;j<nr_class;j++)
2190 int si = start[i], sj = start[j];
2191 int ci = count[i], cj = count[j];
2193 sub_prob.x = Malloc(
svm_node *,sub_prob.l);
2194 sub_prob.y = Malloc(
double,sub_prob.l);
2198 sub_prob.x[k] = x[si+k];
2203 sub_prob.x[ci+k] = x[sj+k];
2204 sub_prob.y[ci+k] = -1;
2207 if(param->probability)
2208 svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);
2210 f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
2212 if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
2213 nonzero[si+k] =
true;
2215 if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
2216 nonzero[sj+k] =
true;
2224 model->nr_class = nr_class;
2226 model->label = Malloc(
int,nr_class);
2227 for(i=0;i<nr_class;i++)
2228 model->label[i] = label[i];
2230 model->rho = Malloc(
double,nr_class*(nr_class-1)/2);
2231 for(i=0;i<nr_class*(nr_class-1)/2;i++)
2232 model->rho[i] = f[i].rho;
2234 if(param->probability)
2236 model->probA = Malloc(
double,nr_class*(nr_class-1)/2);
2237 model->probB = Malloc(
double,nr_class*(nr_class-1)/2);
2238 for(i=0;i<nr_class*(nr_class-1)/2;i++)
2240 model->probA[i] = probA[i];
2241 model->probB[i] = probB[i];
2251 int *nz_count = Malloc(
int,nr_class);
2252 model->nSV = Malloc(
int,nr_class);
2253 for(i=0;i<nr_class;i++)
2256 for(
int j=0;j<count[i];j++)
2257 if(nonzero[start[i]+j])
2262 model->nSV[i] = nSV;
2267 info(
"Total nSV = %d\n",total_sv);
2269 model->l = total_sv;
2270 model->SV = Malloc(
svm_node *,total_sv);
2273 if(nonzero[i]) model->SV[p++] = x[i];
2275 int *nz_start = Malloc(
int,nr_class);
2277 for(i=1;i<nr_class;i++)
2278 nz_start[i] = nz_start[i-1]+nz_count[i-1];
2280 model->sv_coef = Malloc(
double *,nr_class-1);
2281 for(i=0;i<nr_class-1;i++)
2282 model->sv_coef[i] = Malloc(
double,total_sv);
2285 for(i=0;i<nr_class;i++)
2286 for(
int j=i+1;j<nr_class;j++)
2297 int q = nz_start[i];
2301 model->sv_coef[j-1][q++] = f[p].alpha[k];
2305 model->sv_coef[i][q++] = f[p].alpha[ci+k];
2318 for(i=0;i<nr_class*(nr_class-1)/2;i++)
2331 int *fold_start = Malloc(
int,nr_fold+1);
2333 int *perm = Malloc(
int,l);
2338 if((param->svm_type == C_SVC ||
2339 param->svm_type == NU_SVC) && nr_fold < l)
2344 svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
2347 int *fold_count = Malloc(
int,nr_fold);
2349 int *index = Malloc(
int,l);
2352 for (c=0; c<nr_class; c++)
2353 for(i=0;i<count[c];i++)
2355 int j = i+rand()%(count[c]-i);
2356 swap(index[start[c]+j],index[start[c]+i]);
2358 for(i=0;i<nr_fold;i++)
2361 for (c=0; c<nr_class;c++)
2362 fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold;
2365 for (i=1;i<=nr_fold;i++)
2366 fold_start[i] = fold_start[i-1]+fold_count[i-1];
2367 for (c=0; c<nr_class;c++)
2368 for(i=0;i<nr_fold;i++)
2370 int begin = start[c]+i*count[c]/nr_fold;
2371 int end = start[c]+(i+1)*count[c]/nr_fold;
2372 for(
int j=begin;j<end;j++)
2374 perm[fold_start[i]] = index[j];
2379 for (i=1;i<=nr_fold;i++)
2380 fold_start[i] = fold_start[i-1]+fold_count[i-1];
2389 for(i=0;i<l;i++) perm[i]=i;
2392 int j = i+rand()%(l-i);
2393 swap(perm[i],perm[j]);
2395 for(i=0;i<=nr_fold;i++)
2396 fold_start[i]=i*l/nr_fold;
2399 for(i=0;i<nr_fold;i++)
2401 int begin = fold_start[i];
2402 int end = fold_start[i+1];
2406 subprob.l = l-(end-begin);
2407 subprob.x = Malloc(
struct svm_node*,subprob.l);
2408 subprob.y = Malloc(
double,subprob.l);
2411 for(j=0;j<begin;j++)
2413 subprob.x[k] = prob->x[perm[j]];
2414 subprob.y[k] = prob->y[perm[j]];
2419 subprob.x[k] = prob->x[perm[j]];
2420 subprob.y[k] = prob->y[perm[j]];
2423 struct svm_model *submodel = svm_train(&subprob,param);
2424 if(param->probability &&
2425 (param->svm_type == C_SVC || param->svm_type == NU_SVC))
2427 double *prob_estimates=Malloc(
double,svm_get_nr_class(submodel));
2428 for(j=begin;j<end;j++)
2429 target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
2430 free(prob_estimates);
2433 for(j=begin;j<end;j++)
2434 target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
2435 svm_free_and_destroy_model(&submodel);
2444 int svm_get_svm_type(
const svm_model *model)
2446 return model->param.svm_type;
2449 int svm_get_nr_class(
const svm_model *model)
2451 return model->nr_class;
2454 void svm_get_labels(
const svm_model *model,
int* label)
2456 if (model->label != NULL)
2457 for(
int i=0;i<model->nr_class;i++)
2458 label[i] = model->label[i];
2461 double svm_get_svr_probability(
const svm_model *model)
2463 if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
2465 return model->probA[0];
2468 fprintf(stderr,
"Model doesn't contain information for SVR probability inference\n");
2473 double svm_predict_values(
const svm_model *model,
const svm_node *x,
double* dec_values)
2476 if(model->param.svm_type == ONE_CLASS ||
2477 model->param.svm_type == EPSILON_SVR ||
2478 model->param.svm_type == NU_SVR)
2480 double *sv_coef = model->sv_coef[0];
2482 for(i=0;i<model->l;i++)
2483 sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
2484 sum -= model->rho[0];
2487 if(model->param.svm_type == ONE_CLASS)
2488 return (sum>0)?1:-1;
2494 int nr_class = model->nr_class;
2497 double *kvalue = Malloc(
double,l);
2499 kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);
2501 int *start = Malloc(
int,nr_class);
2503 for(i=1;i<nr_class;i++)
2504 start[i] = start[i-1]+model->nSV[i-1];
2506 int *vote = Malloc(
int,nr_class);
2507 for(i=0;i<nr_class;i++)
2511 for(i=0;i<nr_class;i++)
2512 for(
int j=i+1;j<nr_class;j++)
2517 int ci = model->nSV[i];
2518 int cj = model->nSV[j];
2521 double *coef1 = model->sv_coef[j-1];
2522 double *coef2 = model->sv_coef[i];
2524 sum += coef1[si+k] * kvalue[si+k];
2526 sum += coef2[sj+k] * kvalue[sj+k];
2527 sum -= model->rho[p];
2528 dec_values[p] = sum;
2530 if(dec_values[p] > 0)
2537 int vote_max_idx = 0;
2538 for(i=1;i<nr_class;i++)
2539 if(vote[i] > vote[vote_max_idx])
2545 return model->label[vote_max_idx];
2551 int nr_class = model->nr_class;
2553 if(model->param.svm_type == ONE_CLASS ||
2554 model->param.svm_type == EPSILON_SVR ||
2555 model->param.svm_type == NU_SVR)
2556 dec_values = Malloc(
double, 1);
2558 dec_values = Malloc(
double, nr_class*(nr_class-1)/2);
2559 double pred_result = svm_predict_values(model, x, dec_values);
2564 double svm_predict_probability(
2567 if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
2568 model->probA!=NULL && model->probB!=NULL)
2571 int nr_class = model->nr_class;
2572 double *dec_values = Malloc(
double, nr_class*(nr_class-1)/2);
2573 svm_predict_values(model, x, dec_values);
2575 double min_prob=1e-7;
2576 double **pairwise_prob=Malloc(
double *,nr_class);
2577 for(i=0;i<nr_class;i++)
2578 pairwise_prob[i]=Malloc(
double,nr_class);
2580 for(i=0;i<nr_class;i++)
2581 for(
int j=i+1;j<nr_class;j++)
2583 pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
2584 pairwise_prob[j][i]=1-pairwise_prob[i][j];
2587 multiclass_probability(nr_class,pairwise_prob,prob_estimates);
2589 int prob_max_idx = 0;
2590 for(i=1;i<nr_class;i++)
2591 if(prob_estimates[i] > prob_estimates[prob_max_idx])
2593 for(i=0;i<nr_class;i++)
2594 free(pairwise_prob[i]);
2596 free(pairwise_prob);
2597 return model->label[prob_max_idx];
2600 return svm_predict(model, x);
2603 static const char *svm_type_table[] =
2605 "c_svc",
"nu_svc",
"one_class",
"epsilon_svr",
"nu_svr",NULL
2608 static const char *kernel_type_table[]=
2610 "linear",
"polynomial",
"rbf",
"sigmoid",
"precomputed",NULL
2613 int svm_save_model(
const char *model_file_name,
const svm_model *model)
2615 FILE *fp = fopen(model_file_name,
"w");
2616 if(fp==NULL)
return -1;
2618 char *old_locale = strdup(setlocale(LC_ALL, NULL));
2619 setlocale(LC_ALL,
"C");
2623 fprintf(fp,
"svm_type %s\n", svm_type_table[param.svm_type]);
2624 fprintf(fp,
"kernel_type %s\n", kernel_type_table[param.kernel_type]);
2626 if(param.kernel_type == POLY)
2627 fprintf(fp,
"degree %d\n", param.degree);
2629 if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
2630 fprintf(fp,
"gamma %g\n", param.gamma);
2632 if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
2633 fprintf(fp,
"coef0 %g\n", param.coef0);
2635 int nr_class = model->nr_class;
2637 fprintf(fp,
"nr_class %d\n", nr_class);
2638 fprintf(fp,
"total_sv %d\n",l);
2642 for(
int i=0;i<nr_class*(nr_class-1)/2;i++)
2643 fprintf(fp,
" %g",model->rho[i]);
2649 fprintf(fp,
"label");
2650 for(
int i=0;i<nr_class;i++)
2651 fprintf(fp,
" %d",model->label[i]);
2657 fprintf(fp,
"probA");
2658 for(
int i=0;i<nr_class*(nr_class-1)/2;i++)
2659 fprintf(fp,
" %g",model->probA[i]);
2664 fprintf(fp,
"probB");
2665 for(
int i=0;i<nr_class*(nr_class-1)/2;i++)
2666 fprintf(fp,
" %g",model->probB[i]);
2672 fprintf(fp,
"nr_sv");
2673 for(
int i=0;i<nr_class;i++)
2674 fprintf(fp,
" %d",model->nSV[i]);
2678 fprintf(fp,
"SV\n");
2679 const double *
const *sv_coef = model->sv_coef;
2680 const svm_node *
const *SV = model->SV;
2682 for(
int i=0;i<l;i++)
2684 for(
int j=0;j<nr_class-1;j++)
2685 fprintf(fp,
"%.16g ",sv_coef[j][i]);
2689 if(param.kernel_type == PRECOMPUTED)
2690 fprintf(fp,
"0:%d ",(
int)(p->value));
2692 while(p->index != -1)
2694 fprintf(fp,
"%d:%.8g ",p->index,p->value);
2700 setlocale(LC_ALL, old_locale);
2703 if (ferror(fp) != 0 || fclose(fp) != 0)
return -1;
2707 static char *line = NULL;
2708 static int max_line_len;
2710 static char* readline(FILE *input)
2714 if(fgets(line,max_line_len,input) == NULL)
2717 while(strrchr(line,
'\n') == NULL)
2720 line = (
char *) realloc(line,max_line_len);
2721 len = (int) strlen(line);
2722 if(fgets(line+len,max_line_len-len,input) == NULL)
2728 svm_model *svm_load_model(
const char *model_file_name)
2730 FILE *fp = fopen(model_file_name,
"rb");
2731 if(fp==NULL)
return NULL;
2733 char *old_locale = strdup(setlocale(LC_ALL, NULL));
2734 setlocale(LC_ALL,
"C");
2741 model->probA = NULL;
2742 model->probB = NULL;
2743 model->label = NULL;
2749 fscanf(fp,
"%80s",cmd);
2751 if(strcmp(cmd,
"svm_type")==0)
2753 fscanf(fp,
"%80s",cmd);
2755 for(i=0;svm_type_table[i];i++)
2757 if(strcmp(svm_type_table[i],cmd)==0)
2763 if(svm_type_table[i] == NULL)
2765 fprintf(stderr,
"unknown svm type.\n");
2767 setlocale(LC_ALL, old_locale);
2776 else if(strcmp(cmd,
"kernel_type")==0)
2778 fscanf(fp,
"%80s",cmd);
2780 for(i=0;kernel_type_table[i];i++)
2782 if(strcmp(kernel_type_table[i],cmd)==0)
2784 param.kernel_type=i;
2788 if(kernel_type_table[i] == NULL)
2790 fprintf(stderr,
"unknown kernel function.\n");
2792 setlocale(LC_ALL, old_locale);
2801 else if(strcmp(cmd,
"degree")==0)
2802 fscanf(fp,
"%d",¶m.degree);
2803 else if(strcmp(cmd,
"gamma")==0)
2804 fscanf(fp,
"%lf",¶m.gamma);
2805 else if(strcmp(cmd,
"coef0")==0)
2806 fscanf(fp,
"%lf",¶m.coef0);
2807 else if(strcmp(cmd,
"nr_class")==0)
2808 fscanf(fp,
"%d",&model->nr_class);
2809 else if(strcmp(cmd,
"total_sv")==0)
2810 fscanf(fp,
"%d",&model->l);
2811 else if(strcmp(cmd,
"rho")==0)
2813 int n = model->nr_class * (model->nr_class-1)/2;
2814 model->rho = Malloc(
double,n);
2815 for(
int i=0;i<n;i++)
2816 fscanf(fp,
"%lf",&model->rho[i]);
2818 else if(strcmp(cmd,
"label")==0)
2820 int n = model->nr_class;
2821 model->label = Malloc(
int,n);
2822 for(
int i=0;i<n;i++)
2823 fscanf(fp,
"%d",&model->label[i]);
2825 else if(strcmp(cmd,
"probA")==0)
2827 int n = model->nr_class * (model->nr_class-1)/2;
2828 model->probA = Malloc(
double,n);
2829 for(
int i=0;i<n;i++)
2830 fscanf(fp,
"%lf",&model->probA[i]);
2832 else if(strcmp(cmd,
"probB")==0)
2834 int n = model->nr_class * (model->nr_class-1)/2;
2835 model->probB = Malloc(
double,n);
2836 for(
int i=0;i<n;i++)
2837 fscanf(fp,
"%lf",&model->probB[i]);
2839 else if(strcmp(cmd,
"nr_sv")==0)
2841 int n = model->nr_class;
2842 model->nSV = Malloc(
int,n);
2843 for(
int i=0;i<n;i++)
2844 fscanf(fp,
"%d",&model->nSV[i]);
2846 else if(strcmp(cmd,
"SV")==0)
2851 if(c==EOF || c==
'\n')
break;
2857 fprintf(stderr,
"unknown text in model file: [%s]\n",cmd);
2859 setlocale(LC_ALL, old_locale);
2872 long pos = ftell(fp);
2874 max_line_len = 1024;
2875 line = Malloc(
char,max_line_len);
2876 char *p,*endptr,*idx,*val;
2878 while(readline(fp)!=NULL)
2880 p = strtok(line,
":");
2883 p = strtok(NULL,
":");
2889 elements += model->l;
2891 fseek(fp,pos,SEEK_SET);
2893 int m = model->nr_class - 1;
2895 model->sv_coef = Malloc(
double *,m);
2898 model->sv_coef[i] = Malloc(
double,l);
2901 if(l>0) x_space = Malloc(
svm_node,elements);
2907 model->SV[i] = &x_space[j];
2909 p = strtok(line,
" \t");
2910 model->sv_coef[0][i] = strtod(p,&endptr);
2911 for(
int k=1;k<m;k++)
2913 p = strtok(NULL,
" \t");
2914 model->sv_coef[k][i] = strtod(p,&endptr);
2919 idx = strtok(NULL,
":");
2920 val = strtok(NULL,
" \t");
2924 x_space[j].index = (int) strtol(idx,&endptr,10);
2925 x_space[j].value = strtod(val,&endptr);
2929 x_space[j++].index = -1;
2933 setlocale(LC_ALL, old_locale);
2936 if (ferror(fp) != 0 || fclose(fp) != 0)
2943 void svm_free_model_content(
svm_model* model_ptr)
2945 if(model_ptr->free_sv && model_ptr->l > 0 && model_ptr->SV != NULL)
2946 free((
void *)(model_ptr->SV[0]));
2947 if(model_ptr->sv_coef)
2949 for(
int i=0;i<model_ptr->nr_class-1;i++)
2950 free(model_ptr->sv_coef[i]);
2953 free(model_ptr->SV);
2954 model_ptr->SV = NULL;
2956 free(model_ptr->sv_coef);
2957 model_ptr->sv_coef = NULL;
2959 free(model_ptr->rho);
2960 model_ptr->rho = NULL;
2962 free(model_ptr->label);
2963 model_ptr->label= NULL;
2965 free(model_ptr->probA);
2966 model_ptr->probA = NULL;
2968 free(model_ptr->probB);
2969 model_ptr->probB= NULL;
2971 free(model_ptr->nSV);
2972 model_ptr->nSV = NULL;
2975 void svm_free_and_destroy_model(
svm_model** model_ptr_ptr)
2977 if(model_ptr_ptr != NULL && *model_ptr_ptr != NULL)
2979 svm_free_model_content(*model_ptr_ptr);
2980 free(*model_ptr_ptr);
2981 *model_ptr_ptr = NULL;
2987 free(param->weight_label);
2988 free(param->weight);
2995 int svm_type = param->svm_type;
2996 if(svm_type != C_SVC &&
2997 svm_type != NU_SVC &&
2998 svm_type != ONE_CLASS &&
2999 svm_type != EPSILON_SVR &&
3001 return "unknown svm type";
3005 int kernel_type = param->kernel_type;
3006 if(kernel_type != LINEAR &&
3007 kernel_type != POLY &&
3008 kernel_type != RBF &&
3009 kernel_type != SIGMOID &&
3010 kernel_type != PRECOMPUTED)
3011 return "unknown kernel type";
3013 if(param->gamma < 0)
3016 if(param->degree < 0)
3017 return "degree of polynomial kernel < 0";
3021 if(param->cache_size <= 0)
3022 return "cache_size <= 0";
3027 if(svm_type == C_SVC ||
3028 svm_type == EPSILON_SVR ||
3033 if(svm_type == NU_SVC ||
3034 svm_type == ONE_CLASS ||
3036 if(param->nu <= 0 || param->nu > 1)
3037 return "nu <= 0 or nu > 1";
3039 if(svm_type == EPSILON_SVR)
3043 if(param->shrinking != 0 &&
3044 param->shrinking != 1)
3045 return "shrinking != 0 and shrinking != 1";
3047 if(param->probability != 0 &&
3048 param->probability != 1)
3049 return "probability != 0 and probability != 1";
3051 if(param->probability == 1 &&
3052 svm_type == ONE_CLASS)
3053 return "one-class SVM probability output not supported yet";
3058 if(svm_type == NU_SVC)
3061 int max_nr_class = 16;
3063 int *label = Malloc(
int,max_nr_class);
3064 int *count = Malloc(
int,max_nr_class);
3069 int this_label = (int)prob->y[i];
3071 for(j=0;j<nr_class;j++)
3072 if(this_label == label[j])
3079 if(nr_class == max_nr_class)
3082 label = (
int *)realloc(label,max_nr_class*
sizeof(
int));
3083 count = (
int *)realloc(count,max_nr_class*
sizeof(
int));
3085 label[nr_class] = this_label;
3086 count[nr_class] = 1;
3091 for(i=0;i<nr_class;i++)
3094 for(
int j=i+1;j<nr_class;j++)
3097 if(param->nu*(n1+n2)/2 > min(n1,n2))
3101 return "specified nu is infeasible";
3112 int svm_check_probability_model(
const svm_model *model)
3114 return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
3115 model->probA!=NULL && model->probB!=NULL) ||
3116 ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
3117 model->probA!=NULL);
3120 void svm_set_print_string_function(
void (*print_func)(
const char *))
3122 if(print_func == NULL)
3123 svm_print_string = &print_string_stdout;
3125 svm_print_string = print_func;