36 #ifndef DOXYGEN_SHOULD_SKIP_THIS
64 template <
class S,
class T>
inline void clone(T*& dst, S* src, int32_t n)
67 memcpy((
void *)dst,(
void *)src,
sizeof(T)*n);
84 Cache(int32_t l, int64_t size);
90 int32_t get_data(
const int32_t index, Qfloat **data, int32_t len);
91 void swap_index(int32_t i, int32_t j);
105 void lru_delete(head_t *h);
106 void lru_insert(head_t *h);
109 Cache::Cache(int32_t l_, int64_t size_):l(l_),size(size_)
111 head = (head_t *)calloc(l,
sizeof(head_t));
112 size /=
sizeof(Qfloat);
113 size -= l *
sizeof(head_t) /
sizeof(Qfloat);
115 lru_head.next = lru_head.prev = &lru_head;
120 for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
125 void Cache::lru_delete(head_t *h)
128 h->prev->next = h->next;
129 h->next->prev = h->prev;
132 void Cache::lru_insert(head_t *h)
136 h->prev = lru_head.prev;
141 int32_t Cache::get_data(
const int32_t index, Qfloat **data, int32_t len)
143 head_t *h = &head[index];
144 if(h->len) lru_delete(h);
145 int32_t more = len - h->len;
152 head_t *old = lru_head.next;
171 void Cache::swap_index(int32_t i, int32_t j)
175 if(head[i].len) lru_delete(&head[i]);
176 if(head[j].len) lru_delete(&head[j]);
179 if(head[i].len) lru_insert(&head[i]);
180 if(head[j].len) lru_insert(&head[j]);
183 for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
211 virtual Qfloat *get_Q(int32_t column, int32_t len)
const = 0;
212 virtual Qfloat *get_QD()
const = 0;
213 virtual void swap_index(int32_t i, int32_t j)
const = 0;
214 virtual ~QMatrix() {}
222 struct Q_THREAD_PARAM
229 const LibSVMKernel* q;
234 class LibSVMKernel:
public QMatrix {
236 LibSVMKernel(int32_t l, svm_node *
const * x,
const svm_parameter& param);
237 virtual ~LibSVMKernel();
239 virtual Qfloat *get_Q(int32_t column, int32_t len)
const = 0;
240 virtual Qfloat *get_QD()
const = 0;
241 virtual void swap_index(int32_t i, int32_t j)
const
247 static void* compute_Q_parallel_helper(
void* p)
249 Q_THREAD_PARAM* params= (Q_THREAD_PARAM*) p;
251 int32_t start=params->start;
252 int32_t end=params->end;
254 Qfloat* data=params->data;
255 const LibSVMKernel* q=params->q;
259 for(int32_t j=start;j<end;j++)
260 data[j] = (Qfloat) y[i]*y[j]*q->kernel_function(i,j);
264 for(int32_t j=start;j<end;j++)
265 data[j] = (Qfloat) q->kernel_function(i,j);
271 void compute_Q_parallel(Qfloat* data,
float64_t* lab, int32_t i, int32_t start, int32_t len)
const
276 Q_THREAD_PARAM params;
283 compute_Q_parallel_helper((
void*) ¶ms);
287 int32_t total_num=(len-start);
288 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads-1);
289 Q_THREAD_PARAM* params =
SG_MALLOC(Q_THREAD_PARAM, num_threads);
290 int32_t step= total_num/num_threads;
295 for (t=0; t<num_threads; t++)
298 params[t].start=t*step;
299 params[t].end=(t+1)*step;
304 int code=pthread_create(&threads[t], NULL,
305 compute_Q_parallel_helper, (
void*)¶ms[t]);
309 SG_SWARNING(
"Thread creation failed (thread %d of %d) "
310 "with error:'%s'\n",t, num_threads, strerror(code));
317 params[t].start=t*step;
322 compute_Q_parallel_helper(¶ms[t]);
324 for (t=0; t<num_threads; t++)
326 if (pthread_join(threads[t], NULL) != 0)
327 SG_SWARNING(
"pthread_join of thread %d/%d failed\n", t, num_threads);
335 inline float64_t kernel_function(int32_t i, int32_t j)
const
337 return kernel->kernel(x[i]->index,x[j]->index);
346 const int32_t kernel_type;
347 const int32_t degree;
352 LibSVMKernel::LibSVMKernel(int32_t l, svm_node *
const * x_,
const svm_parameter& param)
353 : kernel_type(param.kernel_type), degree(param.degree),
354 gamma(param.gamma), coef0(param.coef0)
362 LibSVMKernel::~LibSVMKernel()
389 virtual ~Solver() {};
391 struct SolutionInfo {
400 int32_t l,
const QMatrix& Q,
const float64_t *p_,
const schar *y_,
402 SolutionInfo* si, int32_t shrinking,
bool use_bias);
408 enum { LOWER_BOUND, UPPER_BOUND, FREE };
423 return (y[i] > 0)? Cp : Cn;
425 void update_alpha_status(int32_t i)
427 if(alpha[i] >= get_C(i))
428 alpha_status[i] = UPPER_BOUND;
429 else if(alpha[i] <= 0)
430 alpha_status[i] = LOWER_BOUND;
431 else alpha_status[i] = FREE;
433 bool is_upper_bound(int32_t i) {
return alpha_status[i] == UPPER_BOUND; }
434 bool is_lower_bound(int32_t i) {
return alpha_status[i] == LOWER_BOUND; }
435 bool is_free(int32_t i) {
return alpha_status[i] == FREE; }
436 void swap_index(int32_t i, int32_t j);
437 void reconstruct_gradient();
438 virtual int32_t select_working_set(int32_t &i, int32_t &j,
float64_t &gap);
440 virtual void do_shrinking();
445 void Solver::swap_index(int32_t i, int32_t j)
457 void Solver::reconstruct_gradient()
461 if(active_size == l)
return;
466 for(j=active_size;j<l;j++)
467 G[j] = G_bar[j] + p[j];
469 for(j=0;j<active_size;j++)
473 if (nr_free*l > 2*active_size*(l-active_size))
475 for(i=active_size;i<l;i++)
477 const Qfloat *Q_i = Q->get_Q(i,active_size);
478 for(j=0;j<active_size;j++)
480 G[i] += alpha[j] * Q_i[j];
485 for(i=0;i<active_size;i++)
488 const Qfloat *Q_i = Q->get_Q(i,l);
490 for(j=active_size;j<l;j++)
491 G[j] += alpha_i * Q_i[j];
497 int32_t p_l,
const QMatrix& p_Q,
const float64_t *p_p,
499 float64_t p_eps, SolutionInfo* p_si, int32_t shrinking,
bool use_bias)
506 clone(alpha,p_alpha,l);
515 for(int32_t i=0;i<l;i++)
516 update_alpha_status(i);
522 for(int32_t i=0;i<l;i++)
539 SG_SINFO(
"Computing gradient for initial set of non-zero alphas\n");
543 if(!is_lower_bound(i))
545 const Qfloat *Q_i = Q->get_Q(i,l);
549 G[j] += alpha_i*Q_i[j];
550 if(is_upper_bound(i))
552 G_bar[j] += get_C(i) * Q_i[j];
566 if (Q->max_train_time > 0 && start_time.cur_time_diff() > Q->max_train_time)
574 if(shrinking) do_shrinking();
580 if(select_working_set(i,j, gap)!=0)
583 reconstruct_gradient();
587 if(select_working_set(i,j, gap)!=0)
599 const Qfloat *Q_i = Q->get_Q(i,active_size);
600 const Qfloat *Q_j = Q->get_Q(j,active_size);
610 double pi=G[i]-Q_i[i]*alpha[i]-Q_i[j]*alpha[j];
611 double pj=G[j]-Q_i[j]*alpha[i]-Q_j[j]*alpha[j];
612 double det=Q_i[i]*Q_j[j]-Q_i[j]*Q_i[j];
613 double alpha_i=-(Q_j[j]*pi-Q_i[j]*pj)/det;
615 double alpha_j=-(-Q_i[j]*pi+Q_i[i]*pj)/det;
618 if (alpha_i==0 || alpha_i == C_i)
620 if (alpha_j==0 || alpha_j == C_j)
623 alpha[i]=alpha_i; alpha[j]=alpha_j;
629 float64_t quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j];
632 float64_t delta = (-G[i]-G[j])/quad_coef;
658 alpha[j] = C_i - diff;
666 alpha[i] = C_j + diff;
672 float64_t quad_coef = Q_i[i]+Q_j[j]-2*Q_i[j];
685 alpha[j] = sum - C_i;
701 alpha[i] = sum - C_j;
717 float64_t delta_alpha_i = alpha[i] - old_alpha_i;
718 float64_t delta_alpha_j = alpha[j] - old_alpha_j;
720 for(int32_t k=0;k<active_size;k++)
722 G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
728 bool ui = is_upper_bound(i);
729 bool uj = is_upper_bound(j);
730 update_alpha_status(i);
731 update_alpha_status(j);
733 if(ui != is_upper_bound(i))
738 G_bar[k] -= C_i * Q_i[k];
741 G_bar[k] += C_i * Q_i[k];
744 if(uj != is_upper_bound(j))
749 G_bar[k] -= C_j * Q_j[k];
752 G_bar[k] += C_j * Q_j[k];
761 v += alpha[i] * (G[i] + p[i]);
767 SG_SPRINT(
"dual obj=%f primal obf=%f gap=%f\n", v/2, primal, gap);
777 p_si->rho = calculate_rho();
784 v += alpha[i] * (G[i] + p[i]);
791 for(int32_t i=0;i<l;i++)
792 p_alpha[active_set[i]] = alpha[i];
795 p_si->upper_bound_p = Cp;
796 p_si->upper_bound_n = Cn;
798 SG_SINFO(
"\noptimization finished, #iter = %d\n",iter);
810 int32_t Solver::select_working_set(
811 int32_t &out_i, int32_t &out_j,
float64_t &gap)
821 int32_t Gmax_idx = -1;
822 int32_t Gmin_idx = -1;
825 for(int32_t t=0;t<active_size;t++)
828 if(!is_upper_bound(t))
837 if(!is_lower_bound(t))
845 int32_t i = Gmax_idx;
846 const Qfloat *Q_i = NULL;
848 Q_i = Q->get_Q(i,active_size);
850 for(int32_t j=0;j<active_size;j++)
854 if (!is_lower_bound(j))
862 float64_t quad_coef=Q_i[i]+QD[j]-2.0*y[i]*Q_i[j];
864 obj_diff = -(grad_diff*grad_diff)/quad_coef;
866 obj_diff = -(grad_diff*grad_diff)/TAU;
868 if (obj_diff <= obj_diff_min)
871 obj_diff_min = obj_diff;
878 if (!is_upper_bound(j))
886 float64_t quad_coef=Q_i[i]+QD[j]+2.0*y[i]*Q_i[j];
888 obj_diff = -(grad_diff*grad_diff)/quad_coef;
890 obj_diff = -(grad_diff*grad_diff)/TAU;
892 if (obj_diff <= obj_diff_min)
895 obj_diff_min = obj_diff;
913 if(is_upper_bound(i))
916 return(-G[i] > Gmax1);
918 return(-G[i] > Gmax2);
920 else if(is_lower_bound(i))
923 return(G[i] > Gmax2);
925 return(G[i] > Gmax1);
932 void Solver::do_shrinking()
939 for(i=0;i<active_size;i++)
943 if(!is_upper_bound(i))
948 if(!is_lower_bound(i))
956 if(!is_upper_bound(i))
961 if(!is_lower_bound(i))
969 if(unshrink ==
false && Gmax1 + Gmax2 <= eps*10)
972 reconstruct_gradient();
976 for(i=0;i<active_size;i++)
977 if (be_shrunk(i, Gmax1, Gmax2))
980 while (active_size > i)
982 if (!be_shrunk(active_size, Gmax1, Gmax2))
984 swap_index(i,active_size);
996 float64_t ub = INF, lb = -INF, sum_free = 0;
997 for(int32_t i=0;i<active_size;i++)
1001 if(is_upper_bound(i))
1008 else if(is_lower_bound(i))
1023 r = sum_free/nr_free;
1034 class WeightedSolver :
public Solver
1042 this->Cs = cost_vec;
1064 class Solver_NU :
public Solver
1069 int32_t p_l,
const QMatrix& p_Q,
const float64_t *p_p,
1071 float64_t p_eps, SolutionInfo* p_si, int32_t shrinking,
bool use_bias)
1074 Solver::Solve(p_l,p_Q,p_p,p_y,p_alpha,p_Cp,p_Cn,p_eps,p_si,
1075 shrinking,use_bias);
1079 int32_t select_working_set(int32_t &i, int32_t &j,
float64_t &gap);
1084 void do_shrinking();
1088 int32_t Solver_NU::select_working_set(
1089 int32_t &out_i, int32_t &out_j,
float64_t &gap)
1099 int32_t Gmaxp_idx = -1;
1103 int32_t Gmaxn_idx = -1;
1105 int32_t Gmin_idx = -1;
1108 for(int32_t t=0;t<active_size;t++)
1111 if(!is_upper_bound(t))
1120 if(!is_lower_bound(t))
1128 int32_t ip = Gmaxp_idx;
1129 int32_t
in = Gmaxn_idx;
1130 const Qfloat *Q_ip = NULL;
1131 const Qfloat *Q_in = NULL;
1133 Q_ip = Q->get_Q(ip,active_size);
1135 Q_in = Q->get_Q(in,active_size);
1137 for(int32_t j=0;j<active_size;j++)
1141 if (!is_lower_bound(j))
1149 float64_t quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];
1151 obj_diff = -(grad_diff*grad_diff)/quad_coef;
1153 obj_diff = -(grad_diff*grad_diff)/TAU;
1155 if (obj_diff <= obj_diff_min)
1158 obj_diff_min = obj_diff;
1165 if (!is_upper_bound(j))
1168 if (-G[j] >= Gmaxn2)
1175 obj_diff = -(grad_diff*grad_diff)/quad_coef;
1177 obj_diff = -(grad_diff*grad_diff)/TAU;
1179 if (obj_diff <= obj_diff_min)
1182 obj_diff_min = obj_diff;
1193 if (y[Gmin_idx] == +1)
1202 bool Solver_NU::be_shrunk(
1206 if(is_upper_bound(i))
1209 return(-G[i] > Gmax1);
1211 return(-G[i] > Gmax4);
1213 else if(is_lower_bound(i))
1216 return(G[i] > Gmax2);
1218 return(G[i] > Gmax3);
1224 void Solver_NU::do_shrinking()
1233 for(i=0;i<active_size;i++)
1235 if(!is_upper_bound(i))
1239 if(-G[i] > Gmax1) Gmax1 = -G[i];
1241 else if(-G[i] > Gmax4) Gmax4 = -G[i];
1243 if(!is_lower_bound(i))
1247 if(G[i] > Gmax2) Gmax2 = G[i];
1249 else if(G[i] > Gmax3) Gmax3 = G[i];
1253 if(unshrink ==
false &&
CMath::max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10)
1256 reconstruct_gradient();
1260 for(i=0;i<active_size;i++)
1261 if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
1264 while (active_size > i)
1266 if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
1268 swap_index(i,active_size);
1278 int32_t nr_free1 = 0,nr_free2 = 0;
1283 for(int32_t i=0; i<active_size; i++)
1287 if(is_upper_bound(i))
1289 else if(is_lower_bound(i))
1299 if(is_upper_bound(i))
1301 else if(is_lower_bound(i))
1313 r1 = sum_free1/nr_free1;
1318 r2 = sum_free2/nr_free2;
1326 class SVC_QMC:
public LibSVMKernel
1329 SVC_QMC(
const svm_problem& prob,
const svm_parameter& param,
const schar *y_, int32_t n_class,
float64_t fac)
1330 :LibSVMKernel(prob.l, prob.x, param)
1335 cache =
new Cache(prob.l,(int64_t)(param.cache_size*(1l<<20)));
1337 for(int32_t i=0;i<prob.l;i++)
1339 QD[i]= factor*(nr_class-1)*kernel_function(i,i);
1343 Qfloat *get_Q(int32_t i, int32_t len)
const
1347 if((start = cache->get_data(i,&data,len)) < len)
1349 compute_Q_parallel(data, NULL, i, start, len);
1351 for(int32_t j=start;j<len;j++)
1354 data[j] *= (factor*(nr_class-1));
1356 data[j] *= (-factor);
1362 inline Qfloat get_orig_Qij(Qfloat Q, int32_t i, int32_t j)
1365 return Q/(nr_class-1);
1370 Qfloat *get_QD()
const
1375 void swap_index(int32_t i, int32_t j)
const
1377 cache->swap_index(i,j);
1378 LibSVMKernel::swap_index(i,j);
1402 class Solver_NUMC :
public Solver
1405 Solver_NUMC(int32_t n_class,
float64_t svm_nu)
1412 int32_t p_l,
const QMatrix& p_Q,
const float64_t *p_p,
1414 float64_t p_eps, SolutionInfo* p_si, int32_t shrinking,
bool use_bias)
1417 Solver::Solve(p_l,p_Q,p_p,p_y,p_alpha,p_Cp,p_Cn,p_eps,p_si,shrinking, use_bias);
1423 int32_t select_working_set(int32_t &i, int32_t &j,
float64_t &gap);
1428 void do_shrinking();
1438 clone(alpha,p_alpha,l);
1441 for(int32_t i=0;i<l;i++)
1442 update_alpha_status(i);
1447 for (int32_t i=0; i<nr_class; i++)
1453 for (int32_t i=0; i<active_size; i++)
1455 update_alpha_status(i);
1456 if(!is_upper_bound(i) && !is_lower_bound(i))
1457 class_count[(int32_t) y[i]]++;
1470 for (int32_t i=0; i<nr_class; i++)
1472 zero_counts[i]=-INF;
1476 for (int32_t i=0; i<active_size; i++)
1482 Qfloat* Q_i = Q->get_Q(i,active_size);
1485 for (
int j=0; j<active_size; j++)
1487 quad+= alpha[i]*alpha[j]*Q_i[j];
1490 if(!is_upper_bound(i) && !is_lower_bound(i))
1495 if (class_count[(int32_t) y[i]] == 0 && y[j]==y[i])
1496 sum_zero_count+= tmp;
1498 SVC_QMC* QMC=(SVC_QMC*) Q;
1499 float64_t norm_tmp=alpha[i]*alpha[j]*QMC->get_orig_Qij(Q_i[j], i, j);
1501 normwcw[(int32_t) y[i]]+=norm_tmp;
1503 normwcw[(int32_t) y[i]]-=2.0/nr_class*norm_tmp;
1504 normwc_const+=norm_tmp;
1507 if (class_count[(int32_t) y[i]] == 0)
1509 if (zero_counts[(int32_t) y[i]]<sum_zero_count)
1510 zero_counts[(int32_t) y[i]]=sum_zero_count;
1513 biases[(int32_t) y[i]+1]-=sum_free;
1514 if (class_count[(int32_t) y[i]] != 0.0)
1515 rho+=sum_free/class_count[(int32_t) y[i]];
1516 outputs[i]+=sum_free+sum_atbound;
1519 for (int32_t i=0; i<nr_class; i++)
1521 if (class_count[i] == 0.0)
1522 rho+=zero_counts[i];
1524 normwcw[i]+=normwc_const/
CMath::sq(nr_class);
1535 for (int32_t i=0; i<nr_class; i++)
1537 if (class_count[i] != 0.0)
1538 biases[i+1]=biases[i+1]/class_count[i]+rho;
1540 biases[i+1]+=rho-zero_counts[i];
1550 for (int32_t i=0; i<l; i++)
1551 outputs[i]+=biases[(int32_t) y[i]+1];
1559 for (int32_t i=0; i<active_size; i++)
1561 if (is_lower_bound(i))
1570 float64_t primal=0.5*quad- nr_class*rho+xi*mu;
1583 int32_t Solver_NUMC::select_working_set(
1584 int32_t &out_i, int32_t &out_j,
float64_t &gap)
1594 int32_t best_out_i=-1;
1595 int32_t best_out_j=-1;
1599 int32_t* Gmaxp_idx =
SG_MALLOC(int32_t, nr_class);
1601 int32_t* Gmin_idx =
SG_MALLOC(int32_t, nr_class);
1604 for (int32_t i=0; i<nr_class; i++)
1610 obj_diff_min[i]=INF;
1613 for(int32_t t=0;t<active_size;t++)
1616 if(!is_upper_bound(t))
1618 if(-G[t] >= Gmaxp[cidx])
1620 Gmaxp[cidx] = -G[t];
1621 Gmaxp_idx[cidx] = t;
1626 for(int32_t j=0;j<active_size;j++)
1629 int32_t ip = Gmaxp_idx[cidx];
1630 const Qfloat *Q_ip = NULL;
1632 Q_ip = Q->get_Q(ip,active_size);
1634 if (!is_lower_bound(j))
1637 if (G[j] >= Gmaxp2[cidx])
1638 Gmaxp2[cidx] = G[j];
1642 float64_t quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];
1644 obj_diff = -(grad_diff*grad_diff)/quad_coef;
1646 obj_diff = -(grad_diff*grad_diff)/TAU;
1648 if (obj_diff <= obj_diff_min[cidx])
1651 obj_diff_min[cidx] = obj_diff;
1656 gap=Gmaxp[cidx]+Gmaxp2[cidx];
1657 if (gap>=best_gap && Gmin_idx[cidx]>=0 &&
1658 Gmaxp_idx[cidx]>=0 && Gmin_idx[cidx]<active_size)
1660 out_i = Gmaxp_idx[cidx];
1661 out_j = Gmin_idx[cidx];
1673 SG_SDEBUG(
"i=%d j=%d best_gap=%f y_i=%f y_j=%f\n", out_i, out_j, gap, y[out_i], y[out_j]);
1688 bool Solver_NUMC::be_shrunk(
1695 void Solver_NUMC::do_shrinking()
1708 class SVC_Q:
public LibSVMKernel
1711 SVC_Q(
const svm_problem& prob,
const svm_parameter& param,
const schar *y_)
1712 :LibSVMKernel(prob.l, prob.x, param)
1715 cache =
new Cache(prob.l,(int64_t)(param.cache_size*(1l<<20)));
1717 for(int32_t i=0;i<prob.l;i++)
1718 QD[i]= (Qfloat)kernel_function(i,i);
1721 Qfloat *get_Q(int32_t i, int32_t len)
const
1725 if((start = cache->get_data(i,&data,len)) < len)
1726 compute_Q_parallel(data, y, i, start, len);
1731 Qfloat *get_QD()
const
1736 void swap_index(int32_t i, int32_t j)
const
1738 cache->swap_index(i,j);
1739 LibSVMKernel::swap_index(i,j);
1757 class ONE_CLASS_Q:
public LibSVMKernel
1760 ONE_CLASS_Q(
const svm_problem& prob,
const svm_parameter& param)
1761 :LibSVMKernel(prob.l, prob.x, param)
1763 cache =
new Cache(prob.l,(int64_t)(param.cache_size*(1l<<20)));
1765 for(int32_t i=0;i<prob.l;i++)
1766 QD[i]= (Qfloat)kernel_function(i,i);
1769 Qfloat *get_Q(int32_t i, int32_t len)
const
1773 if((start = cache->get_data(i,&data,len)) < len)
1774 compute_Q_parallel(data, NULL, i, start, len);
1779 Qfloat *get_QD()
const
1784 void swap_index(int32_t i, int32_t j)
const
1786 cache->swap_index(i,j);
1787 LibSVMKernel::swap_index(i,j);
1801 class SVR_Q:
public LibSVMKernel
1804 SVR_Q(
const svm_problem& prob,
const svm_parameter& param)
1805 :LibSVMKernel(prob.l, prob.x, param)
1808 cache =
new Cache(l,(int64_t)(param.cache_size*(1l<<20)));
1812 for(int32_t k=0;k<l;k++)
1818 QD[k]= (Qfloat)kernel_function(k,k);
1826 void swap_index(int32_t i, int32_t j)
const
1833 Qfloat *get_Q(int32_t i, int32_t len)
const
1836 int32_t real_i = index[i];
1837 if(cache->get_data(real_i,&data,l) < l)
1838 compute_Q_parallel(data, NULL, real_i, 0, l);
1841 Qfloat *buf = buffer[next_buffer];
1842 next_buffer = 1 - next_buffer;
1844 for(int32_t j=0;j<len;j++)
1845 buf[j] = si * sign[j] * data[index[j]];
1849 Qfloat *get_QD()
const
1869 mutable int32_t next_buffer;
1877 static void solve_c_svc(
1878 const svm_problem *prob,
const svm_parameter* param,
1881 int32_t l = prob->l;
1889 if(prob->y[i] > 0) y[i] = +1;
else y[i]=-1;
1893 s.Solve(l, SVC_Q(*prob,*param,y), prob->pv, y,
1894 alpha, Cp, Cn, param->eps, si, param->shrinking, param->use_bias);
1898 sum_alpha += alpha[i];
1901 SG_SINFO(
"nu = %f\n", sum_alpha/(param->C*prob->l));
1911 void solve_c_svc_weighted(
1912 const svm_problem *prob,
const svm_parameter* param,
1925 if(prob->y[i] > 0) y[i] = +1;
else y[i]=-1;
1928 WeightedSolver s = WeightedSolver(prob->C);
1929 s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
1930 alpha, Cp, Cn, param->eps, si, param->shrinking, param->use_bias);
1934 sum_alpha += alpha[i];
1946 static void solve_nu_svc(
1947 const svm_problem *prob,
const svm_parameter *param,
1948 float64_t *alpha, Solver::SolutionInfo* si)
1951 int32_t l = prob->l;
1969 sum_pos -= alpha[i];
1974 sum_neg -= alpha[i];
1983 s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
1984 alpha, 1.0, 1.0, param->eps, si, param->shrinking, param->use_bias);
1994 si->upper_bound_p = 1/r;
1995 si->upper_bound_n = 1/r;
2001 static void solve_nu_multiclass_svc(
const svm_problem *prob,
2002 const svm_parameter *param, Solver::SolutionInfo* si, svm_model* model)
2004 int32_t l = prob->l;
2010 for(int32_t i=0;i<l;i++)
2016 int32_t nr_class=param->nr_class;
2019 for (int32_t j=0; j<nr_class; j++)
2020 sum_class[j] = nu*l/nr_class;
2022 for(int32_t i=0;i<l;i++)
2024 alpha[i] =
CMath::min(1.0,sum_class[int32_t(y[i])]);
2025 sum_class[int32_t(y[i])] -= alpha[i];
2032 for (int32_t i=0;i<l;i++)
2035 Solver_NUMC s(nr_class, nu);
2038 s.Solve(l, Q, zeros, y,
2039 alpha, 1.0, 1.0, param->eps, si, param->shrinking, param->use_bias);
2042 int32_t* class_sv_count=
SG_MALLOC(int32_t, nr_class);
2044 for (int32_t i=0; i<nr_class; i++)
2045 class_sv_count[i]=0;
2047 for (int32_t i=0; i<l; i++)
2050 class_sv_count[(int32_t) y[i]]++;
2056 model->nr_class = nr_class;
2057 model->label = NULL;
2058 model->SV =
SG_MALLOC(svm_node*,nr_class);
2059 model->nSV =
SG_MALLOC(int32_t, nr_class);
2063 for (int32_t i=0; i<nr_class+1; i++)
2066 model->objective = si->obj;
2068 if (param->use_bias)
2070 SG_SDEBUG(
"Computing biases and primal objective\n");
2071 float64_t primal = s.compute_primal(y, alpha, model->rho, model->normwcw);
2072 SG_SINFO(
"Primal = %10.10f\n", primal);
2075 for (int32_t i=0; i<nr_class; i++)
2077 model->nSV[i]=class_sv_count[i];
2078 model->SV[i] =
SG_MALLOC(svm_node,class_sv_count[i]);
2080 class_sv_count[i]=0;
2083 for (int32_t i=0;i<prob->l;i++)
2085 if(fabs(alpha[i]) > 0)
2087 model->SV[(int32_t) y[i]][class_sv_count[(int32_t) y[i]]].index = prob->x[i]->index;
2088 model->sv_coef[(int32_t) y[i]][class_sv_count[(int32_t) y[i]]] = alpha[i];
2089 class_sv_count[(int32_t) y[i]]++;
2098 static void solve_one_class(
2099 const svm_problem *prob,
const svm_parameter *param,
2100 float64_t *alpha, Solver::SolutionInfo* si)
2102 int32_t l = prob->l;
2107 int32_t n = (int32_t)(param->nu*prob->l);
2112 alpha[n] = param->nu * prob->l - n;
2123 s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
2124 alpha, 1.0, 1.0, param->eps, si, param->shrinking, param->use_bias);
2130 static void solve_epsilon_svr(
2131 const svm_problem *prob,
const svm_parameter *param,
2132 float64_t *alpha, Solver::SolutionInfo* si)
2134 int32_t l = prob->l;
2143 linear_term[i] = param->p - prob->y[i];
2147 linear_term[i+l] = param->p + prob->y[i];
2152 s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
2153 alpha2, param->C, param->C, param->eps, si, param->shrinking, param->use_bias);
2158 alpha[i] = alpha2[i] - alpha2[i+l];
2159 sum_alpha += fabs(alpha[i]);
2161 SG_SINFO(
"nu = %f\n",sum_alpha/(param->C*l));
2168 static void solve_nu_svr(
2169 const svm_problem *prob,
const svm_parameter *param,
2170 float64_t *alpha, Solver::SolutionInfo* si)
2172 int32_t l = prob->l;
2185 linear_term[i] = - prob->y[i];
2188 linear_term[i+l] = prob->y[i];
2193 s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
2194 alpha2, C, C, param->eps, si, param->shrinking, param->use_bias);
2199 alpha[i] = alpha2[i] - alpha2[i+l];
2209 struct decision_function
2216 decision_function svm_train_one(
2217 const svm_problem *prob,
const svm_parameter *param,
2221 Solver::SolutionInfo si;
2222 switch(param->svm_type)
2225 solve_c_svc(prob,param,alpha,&si,Cp,Cn);
2228 solve_nu_svc(prob,param,alpha,&si);
2231 solve_one_class(prob,param,alpha,&si);
2234 solve_epsilon_svr(prob,param,alpha,&si);
2237 solve_nu_svr(prob,param,alpha,&si);
2241 SG_SINFO(
"obj = %.16f, rho = %.16f\n",si.obj,si.rho);
2244 if (param->svm_type != ONE_CLASS)
2248 for(int32_t i=0;i<prob->l;i++)
2250 if(fabs(alpha[i]) > 0)
2255 if(fabs(alpha[i]) >= si.upper_bound_p)
2260 if(fabs(alpha[i]) >= si.upper_bound_n)
2265 SG_SINFO(
"nSV = %d, nBSV = %d\n",nSV,nBSV);
2268 decision_function f;
2278 void svm_group_classes(
2279 const svm_problem *prob, int32_t *nr_class_ret, int32_t **label_ret,
2280 int32_t **start_ret, int32_t **count_ret, int32_t *perm)
2282 int32_t l = prob->l;
2283 int32_t max_nr_class = 16;
2284 int32_t nr_class = 0;
2285 int32_t *label =
SG_MALLOC(int32_t, max_nr_class);
2286 int32_t *count =
SG_MALLOC(int32_t, max_nr_class);
2287 int32_t *data_label =
SG_MALLOC(int32_t, l);
2292 int32_t this_label=(int32_t) prob->y[i];
2294 for(j=0;j<nr_class;j++)
2296 if(this_label == label[j])
2305 if(nr_class == max_nr_class)
2308 label=
SG_REALLOC(int32_t, label,max_nr_class);
2309 count=
SG_REALLOC(int32_t, count,max_nr_class);
2311 label[nr_class] = this_label;
2312 count[nr_class] = 1;
2317 int32_t *start =
SG_MALLOC(int32_t, nr_class);
2319 for(i=1;i<nr_class;i++)
2320 start[i] = start[i-1]+count[i-1];
2323 perm[start[data_label[i]]] = i;
2324 ++start[data_label[i]];
2327 for(i=1;i<nr_class;i++)
2328 start[i] = start[i-1]+count[i-1];
2330 *nr_class_ret = nr_class;
2340 svm_model *svm_train(
const svm_problem *prob,
const svm_parameter *param)
2342 svm_model *model =
SG_MALLOC(svm_model,1);
2343 model->param = *param;
2346 if(param->svm_type == ONE_CLASS ||
2347 param->svm_type == EPSILON_SVR ||
2348 param->svm_type == NU_SVR)
2350 SG_SINFO(
"training one class svm or doing epsilon sv regression\n");
2353 model->nr_class = 2;
2354 model->label = NULL;
2357 decision_function f = svm_train_one(prob,param,0,0);
2359 model->rho[0] = f.rho;
2360 model->objective = f.objective;
2364 for(i=0;i<prob->l;i++)
2365 if(fabs(f.alpha[i]) > 0) ++nSV;
2370 for(i=0;i<prob->l;i++)
2371 if(fabs(f.alpha[i]) > 0)
2373 model->SV[j] = prob->x[i];
2374 model->sv_coef[0][j] = f.alpha[i];
2380 else if(param->svm_type == NU_MULTICLASS_SVC)
2382 Solver::SolutionInfo si;
2383 solve_nu_multiclass_svc(prob,param,&si,model);
2384 SG_SINFO(
"obj = %.16f, rho = %.16f\n",si.obj,si.rho);
2389 int32_t l = prob->l;
2391 int32_t *label = NULL;
2392 int32_t *start = NULL;
2393 int32_t *count = NULL;
2397 svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
2405 x[i] = prob->x[perm[i]];
2406 C[i] = prob->C[perm[i]];
2410 pv[i] = prob->pv[perm[i]];
2423 for(i=0;i<nr_class;i++)
2424 weighted_C[i] = param->C;
2425 for(i=0;i<param->nr_weight;i++)
2428 for(j=0;j<nr_class;j++)
2429 if(param->weight_label[i] == label[j])
2432 SG_SWARNING(
"warning: class label %d specified in weight is not found\n", param->weight_label[i]);
2434 weighted_C[j] *= param->weight[i];
2442 decision_function *f =
SG_MALLOC(decision_function,nr_class*(nr_class-1)/2);
2445 for(i=0;i<nr_class;i++)
2446 for(int32_t j=i+1;j<nr_class;j++)
2448 svm_problem sub_prob;
2449 int32_t si = start[i], sj = start[j];
2450 int32_t ci = count[i], cj = count[j];
2452 sub_prob.x =
SG_MALLOC(svm_node *,sub_prob.l);
2460 sub_prob.x[k] = x[si+k];
2462 sub_prob.C[k] = C[si+k];
2463 sub_prob.pv[k] = pv[si+k];
2468 sub_prob.x[ci+k] = x[sj+k];
2469 sub_prob.y[ci+k] = -1;
2470 sub_prob.C[ci+k] = C[sj+k];
2471 sub_prob.pv[ci+k] = pv[sj+k];
2473 sub_prob.y[sub_prob.l]=-1;
2474 sub_prob.C[sub_prob.l]=-1;
2475 sub_prob.pv[sub_prob.l]=-1;
2477 f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
2479 if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
2480 nonzero[si+k] =
true;
2482 if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
2483 nonzero[sj+k] =
true;
2493 model->objective = f[0].objective;
2494 model->nr_class = nr_class;
2496 model->label =
SG_MALLOC(int32_t, nr_class);
2497 for(i=0;i<nr_class;i++)
2498 model->label[i] = label[i];
2501 for(i=0;i<nr_class*(nr_class-1)/2;i++)
2502 model->rho[i] = f[i].rho;
2504 int32_t total_sv = 0;
2505 int32_t *nz_count =
SG_MALLOC(int32_t, nr_class);
2506 model->nSV =
SG_MALLOC(int32_t, nr_class);
2507 for(i=0;i<nr_class;i++)
2510 for(int32_t j=0;j<count[i];j++)
2511 if(nonzero[start[i]+j])
2516 model->nSV[i] = nSV;
2520 SG_SINFO(
"Total nSV = %d\n",total_sv);
2522 model->l = total_sv;
2523 model->SV =
SG_MALLOC(svm_node *,total_sv);
2526 if(nonzero[i]) model->SV[p++] = x[i];
2528 int32_t *nz_start =
SG_MALLOC(int32_t, nr_class);
2530 for(i=1;i<nr_class;i++)
2531 nz_start[i] = nz_start[i-1]+nz_count[i-1];
2534 for(i=0;i<nr_class-1;i++)
2538 for(i=0;i<nr_class;i++)
2539 for(int32_t j=i+1;j<nr_class;j++)
2545 int32_t si = start[i];
2546 int32_t sj = start[j];
2547 int32_t ci = count[i];
2548 int32_t cj = count[j];
2550 int32_t q = nz_start[i];
2554 model->sv_coef[j-1][q++] = f[p].alpha[k];
2558 model->sv_coef[i][q++] = f[p].alpha[ci+k];
2571 for(i=0;i<nr_class*(nr_class-1)/2;i++)
2580 void svm_destroy_model(svm_model* model)
2582 if(model->free_sv && model->l > 0)
2583 SG_FREE((
void *)(model->SV[0]));
2584 for(int32_t i=0;i<model->nr_class-1;i++)
2594 void svm_destroy_param(svm_parameter* param)
2600 const char *svm_check_parameter(
2601 const svm_problem *prob,
const svm_parameter *param)
2605 int32_t svm_type = param->svm_type;
2606 if(svm_type != C_SVC &&
2607 svm_type != NU_SVC &&
2608 svm_type != ONE_CLASS &&
2609 svm_type != EPSILON_SVR &&
2610 svm_type != NU_SVR &&
2611 svm_type != NU_MULTICLASS_SVC)
2612 return "unknown svm type";
2616 int32_t kernel_type = param->kernel_type;
2617 if(kernel_type != LINEAR &&
2618 kernel_type != POLY &&
2619 kernel_type != RBF &&
2620 kernel_type != SIGMOID &&
2621 kernel_type != PRECOMPUTED)
2622 return "unknown kernel type";
2624 if(param->degree < 0)
2625 return "degree of polynomial kernel < 0";
2629 if(param->cache_size <= 0)
2630 return "cache_size <= 0";
2635 if(svm_type == C_SVC ||
2636 svm_type == EPSILON_SVR ||
2641 if(svm_type == NU_SVC ||
2642 svm_type == ONE_CLASS ||
2644 if(param->nu <= 0 || param->nu > 1)
2645 return "nu <= 0 or nu > 1";
2647 if(svm_type == EPSILON_SVR)
2651 if(param->shrinking != 0 &&
2652 param->shrinking != 1)
2653 return "shrinking != 0 and shrinking != 1";
2658 if(svm_type == NU_SVC)
2660 int32_t l = prob->l;
2661 int32_t max_nr_class = 16;
2662 int32_t nr_class = 0;
2663 int32_t *label =
SG_MALLOC(int32_t, max_nr_class);
2664 int32_t *count =
SG_MALLOC(int32_t, max_nr_class);
2669 int32_t this_label = (int32_t) prob->y[i];
2671 for(j=0;j<nr_class;j++)
2672 if(this_label == label[j])
2679 if(nr_class == max_nr_class)
2682 label=
SG_REALLOC(int32_t, label, max_nr_class);
2683 count=
SG_REALLOC(int32_t, count, max_nr_class);
2685 label[nr_class] = this_label;
2686 count[nr_class] = 1;
2691 for(i=0;i<nr_class;i++)
2693 int32_t n1 = count[i];
2694 for(int32_t j=i+1;j<nr_class;j++)
2696 int32_t n2 = count[j];
2701 return "specified nu is infeasible";
2712 #endif // DOXYGEN_SHOULD_SKIP_THIS