00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00037
00038 #include "lib/memory.h"
00039 #include "classifier/svm/SVM_libsvm.h"
00040 #include "kernel/Kernel.h"
00041 #include "lib/io.h"
00042 #include "lib/common.h"
00043 #include "lib/Mathematics.h"
00044
00045 #include <stdio.h>
00046 #include <stdlib.h>
00047 #include <ctype.h>
00048 #include <float.h>
00049 #include <string.h>
00050 #include <stdarg.h>
00051 typedef KERNELCACHE_ELEM Qfloat;
00052 typedef float64_t schar;
00053 #ifndef min
00054 template <class T> inline T min(T x,T y) { return (x<y)?x:y; }
00055 #endif
00056 #ifndef max
00057 template <class T> inline T max(T x,T y) { return (x>y)?x:y; }
00058 #endif
00059 template <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
00060 template <class S, class T> inline void clone(T*& dst, S* src, int32_t n)
00061 {
00062 dst = new T[n];
00063 memcpy((void *)dst,(void *)src,sizeof(T)*n);
00064 }
00065 #define INF HUGE_VAL
00066 #define TAU 1e-12
00067 #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
00068
00069
00070
00071
00072
00073
00074
00075 class Cache
00076 {
00077 public:
00078 Cache(int32_t l, int64_t size);
00079 ~Cache();
00080
00081
00082
00083
00084 int32_t get_data(const int32_t index, Qfloat **data, int32_t len);
00085 void swap_index(int32_t i, int32_t j);
00086
00087 private:
00088 int32_t l;
00089 int64_t size;
00090 struct head_t
00091 {
00092 head_t *prev, *next;
00093 Qfloat *data;
00094 int32_t len;
00095 };
00096
00097 head_t *head;
00098 head_t lru_head;
00099 void lru_delete(head_t *h);
00100 void lru_insert(head_t *h);
00101 };
00102
00103 Cache::Cache(int32_t l_, int64_t size_):l(l_),size(size_)
00104 {
00105 head = (head_t *)calloc(l,sizeof(head_t));
00106 size /= sizeof(Qfloat);
00107 size -= l * sizeof(head_t) / sizeof(Qfloat);
00108 size = max(size, (int64_t) 2*l);
00109 lru_head.next = lru_head.prev = &lru_head;
00110 }
00111
00112 Cache::~Cache()
00113 {
00114 for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
00115 free(h->data);
00116 free(head);
00117 }
00118
00119 void Cache::lru_delete(head_t *h)
00120 {
00121
00122 h->prev->next = h->next;
00123 h->next->prev = h->prev;
00124 }
00125
00126 void Cache::lru_insert(head_t *h)
00127 {
00128
00129 h->next = &lru_head;
00130 h->prev = lru_head.prev;
00131 h->prev->next = h;
00132 h->next->prev = h;
00133 }
00134
00135 int32_t Cache::get_data(const int32_t index, Qfloat **data, int32_t len)
00136 {
00137 head_t *h = &head[index];
00138 if(h->len) lru_delete(h);
00139 int32_t more = len - h->len;
00140
00141 if(more > 0)
00142 {
00143
00144 while(size < more)
00145 {
00146 head_t *old = lru_head.next;
00147 lru_delete(old);
00148 free(old->data);
00149 size += old->len;
00150 old->data = 0;
00151 old->len = 0;
00152 }
00153
00154
00155 h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
00156 size -= more;
00157 swap(h->len,len);
00158 }
00159
00160 lru_insert(h);
00161 *data = h->data;
00162 return len;
00163 }
00164
00165 void Cache::swap_index(int32_t i, int32_t j)
00166 {
00167 if(i==j) return;
00168
00169 if(head[i].len) lru_delete(&head[i]);
00170 if(head[j].len) lru_delete(&head[j]);
00171 swap(head[i].data,head[j].data);
00172 swap(head[i].len,head[j].len);
00173 if(head[i].len) lru_insert(&head[i]);
00174 if(head[j].len) lru_insert(&head[j]);
00175
00176 if(i>j) swap(i,j);
00177 for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
00178 {
00179 if(h->len > i)
00180 {
00181 if(h->len > j)
00182 swap(h->data[i],h->data[j]);
00183 else
00184 {
00185
00186 lru_delete(h);
00187 free(h->data);
00188 size += h->len;
00189 h->data = 0;
00190 h->len = 0;
00191 }
00192 }
00193 }
00194 }
00195
00196
00197
00198
00199
00200
00201
00202
00203 class QMatrix {
00204 public:
00205 virtual Qfloat *get_Q(int32_t column, int32_t len) const = 0;
00206 virtual Qfloat *get_QD() const = 0;
00207 virtual void swap_index(int32_t i, int32_t j) const = 0;
00208 virtual ~QMatrix() {}
00209 };
00210
00211 class LibSVMKernel: public QMatrix {
00212 public:
00213 LibSVMKernel(int32_t l, svm_node * const * x, const svm_parameter& param);
00214 virtual ~LibSVMKernel();
00215
00216 virtual Qfloat *get_Q(int32_t column, int32_t len) const = 0;
00217 virtual Qfloat *get_QD() const = 0;
00218 virtual void swap_index(int32_t i, int32_t j) const
00219 {
00220 swap(x[i],x[j]);
00221 if(x_square) swap(x_square[i],x_square[j]);
00222 }
00223
00224 inline float64_t kernel_function(int32_t i, int32_t j) const
00225 {
00226 return kernel->kernel(x[i]->index,x[j]->index);
00227 }
00228
00229 private:
00230 CKernel* kernel;
00231 const svm_node **x;
00232 float64_t *x_square;
00233
00234
00235 const int32_t kernel_type;
00236 const int32_t degree;
00237 const float64_t gamma;
00238 const float64_t coef0;
00239 };
00240
00241 LibSVMKernel::LibSVMKernel(int32_t l, svm_node * const * x_, const svm_parameter& param)
00242 : kernel_type(param.kernel_type), degree(param.degree),
00243 gamma(param.gamma), coef0(param.coef0)
00244 {
00245 clone(x,x_,l);
00246 x_square = 0;
00247 kernel=param.kernel;
00248 }
00249
00250 LibSVMKernel::~LibSVMKernel()
00251 {
00252 delete[] x;
00253 delete[] x_square;
00254 }
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274 class Solver {
00275 public:
00276 Solver() {};
00277 virtual ~Solver() {};
00278
00279 struct SolutionInfo {
00280 float64_t obj;
00281 float64_t rho;
00282 float64_t upper_bound_p;
00283 float64_t upper_bound_n;
00284 float64_t r;
00285 };
00286
00287 void Solve(
00288 int32_t l, const QMatrix& Q, const float64_t *p_, const schar *y_,
00289 float64_t *alpha_, float64_t Cp, float64_t Cn, float64_t eps,
00290 SolutionInfo* si, int32_t shrinking);
00291
00292 protected:
00293 int32_t active_size;
00294 schar *y;
00295 float64_t *G;
00296 enum { LOWER_BOUND, UPPER_BOUND, FREE };
00297 char *alpha_status;
00298 float64_t *alpha;
00299 const QMatrix *Q;
00300 const Qfloat *QD;
00301 float64_t eps;
00302 float64_t Cp,Cn;
00303 float64_t *p;
00304 int32_t *active_set;
00305 float64_t *G_bar;
00306 int32_t l;
00307 bool unshrink;
00308
00309 float64_t get_C(int32_t i)
00310 {
00311 return (y[i] > 0)? Cp : Cn;
00312 }
00313 void update_alpha_status(int32_t i)
00314 {
00315 if(alpha[i] >= get_C(i))
00316 alpha_status[i] = UPPER_BOUND;
00317 else if(alpha[i] <= 0)
00318 alpha_status[i] = LOWER_BOUND;
00319 else alpha_status[i] = FREE;
00320 }
00321 bool is_upper_bound(int32_t i) { return alpha_status[i] == UPPER_BOUND; }
00322 bool is_lower_bound(int32_t i) { return alpha_status[i] == LOWER_BOUND; }
00323 bool is_free(int32_t i) { return alpha_status[i] == FREE; }
00324 void swap_index(int32_t i, int32_t j);
00325 void reconstruct_gradient();
00326 virtual int32_t select_working_set(int32_t &i, int32_t &j, float64_t &gap);
00327 virtual float64_t calculate_rho();
00328 virtual void do_shrinking();
00329 private:
00330 bool be_shrunk(int32_t i, float64_t Gmax1, float64_t Gmax2);
00331 };
00332
00333 void Solver::swap_index(int32_t i, int32_t j)
00334 {
00335 Q->swap_index(i,j);
00336 swap(y[i],y[j]);
00337 swap(G[i],G[j]);
00338 swap(alpha_status[i],alpha_status[j]);
00339 swap(alpha[i],alpha[j]);
00340 swap(p[i],p[j]);
00341 swap(active_set[i],active_set[j]);
00342 swap(G_bar[i],G_bar[j]);
00343 }
00344
00345 void Solver::reconstruct_gradient()
00346 {
00347
00348
00349 if(active_size == l) return;
00350
00351 int32_t i,j;
00352 int32_t nr_free = 0;
00353
00354 for(j=active_size;j<l;j++)
00355 G[j] = G_bar[j] + p[j];
00356
00357 for(j=0;j<active_size;j++)
00358 if(is_free(j))
00359 nr_free++;
00360
00361 if (nr_free*l > 2*active_size*(l-active_size))
00362 {
00363 for(i=active_size;i<l;i++)
00364 {
00365 const Qfloat *Q_i = Q->get_Q(i,active_size);
00366 for(j=0;j<active_size;j++)
00367 if(is_free(j))
00368 G[i] += alpha[j] * Q_i[j];
00369 }
00370 }
00371 else
00372 {
00373 for(i=0;i<active_size;i++)
00374 if(is_free(i))
00375 {
00376 const Qfloat *Q_i = Q->get_Q(i,l);
00377 float64_t alpha_i = alpha[i];
00378 for(j=active_size;j<l;j++)
00379 G[j] += alpha_i * Q_i[j];
00380 }
00381 }
00382 }
00383
00384 void Solver::Solve(
00385 int32_t p_l, const QMatrix& p_Q, const float64_t *p_p,
00386 const schar *p_y, float64_t *p_alpha, float64_t p_Cp, float64_t p_Cn,
00387 float64_t p_eps, SolutionInfo* p_si, int32_t shrinking)
00388 {
00389 this->l = p_l;
00390 this->Q = &p_Q;
00391 QD=Q->get_QD();
00392 clone(p, p_p,l);
00393 clone(y, p_y,l);
00394 clone(alpha,p_alpha,l);
00395 this->Cp = p_Cp;
00396 this->Cn = p_Cn;
00397 this->eps = p_eps;
00398 unshrink = false;
00399
00400
00401 {
00402 alpha_status = new char[l];
00403 for(int32_t i=0;i<l;i++)
00404 update_alpha_status(i);
00405 }
00406
00407
00408 {
00409 active_set = new int32_t[l];
00410 for(int32_t i=0;i<l;i++)
00411 active_set[i] = i;
00412 active_size = l;
00413 }
00414
00415
00416 {
00417 G = new float64_t[l];
00418 G_bar = new float64_t[l];
00419 int32_t i;
00420 for(i=0;i<l;i++)
00421 {
00422 G[i] = p_p[i];
00423 G_bar[i] = 0;
00424 }
00425 for(i=0;i<l;i++)
00426 if(!is_lower_bound(i))
00427 {
00428 const Qfloat *Q_i = Q->get_Q(i,l);
00429 float64_t alpha_i = alpha[i];
00430 int32_t j;
00431 for(j=0;j<l;j++)
00432 G[j] += alpha_i*Q_i[j];
00433 if(is_upper_bound(i))
00434 for(j=0;j<l;j++)
00435 G_bar[j] += get_C(i) * Q_i[j];
00436 }
00437 }
00438
00439
00440
00441 int32_t iter = 0;
00442 int32_t counter = min(l,1000)+1;
00443
00444 while(1)
00445 {
00446
00447
00448 if(--counter == 0)
00449 {
00450 counter = min(l,1000);
00451 if(shrinking) do_shrinking();
00452
00453 }
00454
00455 int32_t i,j;
00456 float64_t gap;
00457 if(select_working_set(i,j, gap)!=0)
00458 {
00459
00460 reconstruct_gradient();
00461
00462 active_size = l;
00463
00464 if(select_working_set(i,j, gap)!=0)
00465 break;
00466 else
00467 counter = 1;
00468 }
00469
00470 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(eps), 6);
00471
00472 ++iter;
00473
00474
00475
00476 const Qfloat *Q_i = Q->get_Q(i,active_size);
00477 const Qfloat *Q_j = Q->get_Q(j,active_size);
00478
00479 float64_t C_i = get_C(i);
00480 float64_t C_j = get_C(j);
00481
00482 float64_t old_alpha_i = alpha[i];
00483 float64_t old_alpha_j = alpha[j];
00484
00485 if(y[i]!=y[j])
00486 {
00487 float64_t quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j];
00488 if (quad_coef <= 0)
00489 quad_coef = TAU;
00490 float64_t delta = (-G[i]-G[j])/quad_coef;
00491 float64_t diff = alpha[i] - alpha[j];
00492 alpha[i] += delta;
00493 alpha[j] += delta;
00494
00495 if(diff > 0)
00496 {
00497 if(alpha[j] < 0)
00498 {
00499 alpha[j] = 0;
00500 alpha[i] = diff;
00501 }
00502 }
00503 else
00504 {
00505 if(alpha[i] < 0)
00506 {
00507 alpha[i] = 0;
00508 alpha[j] = -diff;
00509 }
00510 }
00511 if(diff > C_i - C_j)
00512 {
00513 if(alpha[i] > C_i)
00514 {
00515 alpha[i] = C_i;
00516 alpha[j] = C_i - diff;
00517 }
00518 }
00519 else
00520 {
00521 if(alpha[j] > C_j)
00522 {
00523 alpha[j] = C_j;
00524 alpha[i] = C_j + diff;
00525 }
00526 }
00527 }
00528 else
00529 {
00530 float64_t quad_coef = Q_i[i]+Q_j[j]-2*Q_i[j];
00531 if (quad_coef <= 0)
00532 quad_coef = TAU;
00533 float64_t delta = (G[i]-G[j])/quad_coef;
00534 float64_t sum = alpha[i] + alpha[j];
00535 alpha[i] -= delta;
00536 alpha[j] += delta;
00537
00538 if(sum > C_i)
00539 {
00540 if(alpha[i] > C_i)
00541 {
00542 alpha[i] = C_i;
00543 alpha[j] = sum - C_i;
00544 }
00545 }
00546 else
00547 {
00548 if(alpha[j] < 0)
00549 {
00550 alpha[j] = 0;
00551 alpha[i] = sum;
00552 }
00553 }
00554 if(sum > C_j)
00555 {
00556 if(alpha[j] > C_j)
00557 {
00558 alpha[j] = C_j;
00559 alpha[i] = sum - C_j;
00560 }
00561 }
00562 else
00563 {
00564 if(alpha[i] < 0)
00565 {
00566 alpha[i] = 0;
00567 alpha[j] = sum;
00568 }
00569 }
00570 }
00571
00572
00573
00574 float64_t delta_alpha_i = alpha[i] - old_alpha_i;
00575 float64_t delta_alpha_j = alpha[j] - old_alpha_j;
00576
00577 for(int32_t k=0;k<active_size;k++)
00578 {
00579 G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
00580 }
00581
00582
00583
00584 {
00585 bool ui = is_upper_bound(i);
00586 bool uj = is_upper_bound(j);
00587 update_alpha_status(i);
00588 update_alpha_status(j);
00589 int32_t k;
00590 if(ui != is_upper_bound(i))
00591 {
00592 Q_i = Q->get_Q(i,l);
00593 if(ui)
00594 for(k=0;k<l;k++)
00595 G_bar[k] -= C_i * Q_i[k];
00596 else
00597 for(k=0;k<l;k++)
00598 G_bar[k] += C_i * Q_i[k];
00599 }
00600
00601 if(uj != is_upper_bound(j))
00602 {
00603 Q_j = Q->get_Q(j,l);
00604 if(uj)
00605 for(k=0;k<l;k++)
00606 G_bar[k] -= C_j * Q_j[k];
00607 else
00608 for(k=0;k<l;k++)
00609 G_bar[k] += C_j * Q_j[k];
00610 }
00611 }
00612 }
00613
00614
00615
00616 p_si->rho = calculate_rho();
00617
00618
00619 {
00620 float64_t v = 0;
00621 int32_t i;
00622 for(i=0;i<l;i++)
00623 v += alpha[i] * (G[i] + p[i]);
00624
00625 p_si->obj = v/2;
00626 }
00627
00628
00629 {
00630 for(int32_t i=0;i<l;i++)
00631 p_alpha[active_set[i]] = alpha[i];
00632 }
00633
00634 p_si->upper_bound_p = Cp;
00635 p_si->upper_bound_n = Cn;
00636
00637 SG_SINFO("\noptimization finished, #iter = %d\n",iter);
00638
00639 delete[] p;
00640 delete[] y;
00641 delete[] alpha;
00642 delete[] alpha_status;
00643 delete[] active_set;
00644 delete[] G;
00645 delete[] G_bar;
00646 }
00647
00648
00649 int32_t Solver::select_working_set(
00650 int32_t &out_i, int32_t &out_j, float64_t &gap)
00651 {
00652
00653
00654
00655
00656
00657
00658 float64_t Gmax = -INF;
00659 float64_t Gmax2 = -INF;
00660 int32_t Gmax_idx = -1;
00661 int32_t Gmin_idx = -1;
00662 float64_t obj_diff_min = INF;
00663
00664 for(int32_t t=0;t<active_size;t++)
00665 if(y[t]==+1)
00666 {
00667 if(!is_upper_bound(t))
00668 if(-G[t] >= Gmax)
00669 {
00670 Gmax = -G[t];
00671 Gmax_idx = t;
00672 }
00673 }
00674 else
00675 {
00676 if(!is_lower_bound(t))
00677 if(G[t] >= Gmax)
00678 {
00679 Gmax = G[t];
00680 Gmax_idx = t;
00681 }
00682 }
00683
00684 int32_t i = Gmax_idx;
00685 const Qfloat *Q_i = NULL;
00686 if(i != -1)
00687 Q_i = Q->get_Q(i,active_size);
00688
00689 for(int32_t j=0;j<active_size;j++)
00690 {
00691 if(y[j]==+1)
00692 {
00693 if (!is_lower_bound(j))
00694 {
00695 float64_t grad_diff=Gmax+G[j];
00696 if (G[j] >= Gmax2)
00697 Gmax2 = G[j];
00698 if (grad_diff > 0)
00699 {
00700 float64_t obj_diff;
00701 float64_t quad_coef=Q_i[i]+QD[j]-2.0*y[i]*Q_i[j];
00702 if (quad_coef > 0)
00703 obj_diff = -(grad_diff*grad_diff)/quad_coef;
00704 else
00705 obj_diff = -(grad_diff*grad_diff)/TAU;
00706
00707 if (obj_diff <= obj_diff_min)
00708 {
00709 Gmin_idx=j;
00710 obj_diff_min = obj_diff;
00711 }
00712 }
00713 }
00714 }
00715 else
00716 {
00717 if (!is_upper_bound(j))
00718 {
00719 float64_t grad_diff= Gmax-G[j];
00720 if (-G[j] >= Gmax2)
00721 Gmax2 = -G[j];
00722 if (grad_diff > 0)
00723 {
00724 float64_t obj_diff;
00725 float64_t quad_coef=Q_i[i]+QD[j]+2.0*y[i]*Q_i[j];
00726 if (quad_coef > 0)
00727 obj_diff = -(grad_diff*grad_diff)/quad_coef;
00728 else
00729 obj_diff = -(grad_diff*grad_diff)/TAU;
00730
00731 if (obj_diff <= obj_diff_min)
00732 {
00733 Gmin_idx=j;
00734 obj_diff_min = obj_diff;
00735 }
00736 }
00737 }
00738 }
00739 }
00740
00741 gap=Gmax+Gmax2;
00742 if(gap < eps)
00743 return 1;
00744
00745 out_i = Gmax_idx;
00746 out_j = Gmin_idx;
00747 return 0;
00748 }
00749
00750 bool Solver::be_shrunk(int32_t i, float64_t Gmax1, float64_t Gmax2)
00751 {
00752 if(is_upper_bound(i))
00753 {
00754 if(y[i]==+1)
00755 return(-G[i] > Gmax1);
00756 else
00757 return(-G[i] > Gmax2);
00758 }
00759 else if(is_lower_bound(i))
00760 {
00761 if(y[i]==+1)
00762 return(G[i] > Gmax2);
00763 else
00764 return(G[i] > Gmax1);
00765 }
00766 else
00767 return(false);
00768 }
00769
00770
00771 void Solver::do_shrinking()
00772 {
00773 int32_t i;
00774 float64_t Gmax1 = -INF;
00775 float64_t Gmax2 = -INF;
00776
00777
00778 for(i=0;i<active_size;i++)
00779 {
00780 if(y[i]==+1)
00781 {
00782 if(!is_upper_bound(i))
00783 {
00784 if(-G[i] >= Gmax1)
00785 Gmax1 = -G[i];
00786 }
00787 if(!is_lower_bound(i))
00788 {
00789 if(G[i] >= Gmax2)
00790 Gmax2 = G[i];
00791 }
00792 }
00793 else
00794 {
00795 if(!is_upper_bound(i))
00796 {
00797 if(-G[i] >= Gmax2)
00798 Gmax2 = -G[i];
00799 }
00800 if(!is_lower_bound(i))
00801 {
00802 if(G[i] >= Gmax1)
00803 Gmax1 = G[i];
00804 }
00805 }
00806 }
00807
00808 if(unshrink == false && Gmax1 + Gmax2 <= eps*10)
00809 {
00810 unshrink = true;
00811 reconstruct_gradient();
00812 active_size = l;
00813 }
00814
00815 for(i=0;i<active_size;i++)
00816 if (be_shrunk(i, Gmax1, Gmax2))
00817 {
00818 active_size--;
00819 while (active_size > i)
00820 {
00821 if (!be_shrunk(active_size, Gmax1, Gmax2))
00822 {
00823 swap_index(i,active_size);
00824 break;
00825 }
00826 active_size--;
00827 }
00828 }
00829 }
00830
00831 float64_t Solver::calculate_rho()
00832 {
00833 float64_t r;
00834 int32_t nr_free = 0;
00835 float64_t ub = INF, lb = -INF, sum_free = 0;
00836 for(int32_t i=0;i<active_size;i++)
00837 {
00838 float64_t yG = y[i]*G[i];
00839
00840 if(is_upper_bound(i))
00841 {
00842 if(y[i]==-1)
00843 ub = min(ub,yG);
00844 else
00845 lb = max(lb,yG);
00846 }
00847 else if(is_lower_bound(i))
00848 {
00849 if(y[i]==+1)
00850 ub = min(ub,yG);
00851 else
00852 lb = max(lb,yG);
00853 }
00854 else
00855 {
00856 ++nr_free;
00857 sum_free += yG;
00858 }
00859 }
00860
00861 if(nr_free>0)
00862 r = sum_free/nr_free;
00863 else
00864 r = (ub+lb)/2;
00865
00866 return r;
00867 }
00868
00869
00870
00871
00872
00873
00874 class Solver_NU : public Solver
00875 {
00876 public:
00877 Solver_NU() {}
00878 void Solve(
00879 int32_t p_l, const QMatrix& p_Q, const float64_t *p_p,
00880 const schar *p_y, float64_t* p_alpha, float64_t p_Cp, float64_t p_Cn,
00881 float64_t p_eps, SolutionInfo* p_si, int32_t shrinking)
00882 {
00883 this->si = p_si;
00884 Solver::Solve(p_l,p_Q,p,p_y,p_alpha,p_Cp,p_Cn,p_eps,p_si,shrinking);
00885 }
00886 private:
00887 SolutionInfo *si;
00888 int32_t select_working_set(int32_t &i, int32_t &j, float64_t &gap);
00889 float64_t calculate_rho();
00890 bool be_shrunk(
00891 int32_t i, float64_t Gmax1, float64_t Gmax2, float64_t Gmax3,
00892 float64_t Gmax4);
00893 void do_shrinking();
00894 };
00895
00896
00897 int32_t Solver_NU::select_working_set(
00898 int32_t &out_i, int32_t &out_j, float64_t &gap)
00899 {
00900
00901
00902
00903
00904
00905
00906 float64_t Gmaxp = -INF;
00907 float64_t Gmaxp2 = -INF;
00908 int32_t Gmaxp_idx = -1;
00909
00910 float64_t Gmaxn = -INF;
00911 float64_t Gmaxn2 = -INF;
00912 int32_t Gmaxn_idx = -1;
00913
00914 int32_t Gmin_idx = -1;
00915 float64_t obj_diff_min = INF;
00916
00917 for(int32_t t=0;t<active_size;t++)
00918 if(y[t]==+1)
00919 {
00920 if(!is_upper_bound(t))
00921 if(-G[t] >= Gmaxp)
00922 {
00923 Gmaxp = -G[t];
00924 Gmaxp_idx = t;
00925 }
00926 }
00927 else
00928 {
00929 if(!is_lower_bound(t))
00930 if(G[t] >= Gmaxn)
00931 {
00932 Gmaxn = G[t];
00933 Gmaxn_idx = t;
00934 }
00935 }
00936
00937 int32_t ip = Gmaxp_idx;
00938 int32_t in = Gmaxn_idx;
00939 const Qfloat *Q_ip = NULL;
00940 const Qfloat *Q_in = NULL;
00941 if(ip != -1)
00942 Q_ip = Q->get_Q(ip,active_size);
00943 if(in != -1)
00944 Q_in = Q->get_Q(in,active_size);
00945
00946 for(int32_t j=0;j<active_size;j++)
00947 {
00948 if(y[j]==+1)
00949 {
00950 if (!is_lower_bound(j))
00951 {
00952 float64_t grad_diff=Gmaxp+G[j];
00953 if (G[j] >= Gmaxp2)
00954 Gmaxp2 = G[j];
00955 if (grad_diff > 0)
00956 {
00957 float64_t obj_diff;
00958 float64_t quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];
00959 if (quad_coef > 0)
00960 obj_diff = -(grad_diff*grad_diff)/quad_coef;
00961 else
00962 obj_diff = -(grad_diff*grad_diff)/TAU;
00963
00964 if (obj_diff <= obj_diff_min)
00965 {
00966 Gmin_idx=j;
00967 obj_diff_min = obj_diff;
00968 }
00969 }
00970 }
00971 }
00972 else
00973 {
00974 if (!is_upper_bound(j))
00975 {
00976 float64_t grad_diff=Gmaxn-G[j];
00977 if (-G[j] >= Gmaxn2)
00978 Gmaxn2 = -G[j];
00979 if (grad_diff > 0)
00980 {
00981 float64_t obj_diff;
00982 float64_t quad_coef = Q_in[in]+QD[j]-2*Q_in[j];
00983 if (quad_coef > 0)
00984 obj_diff = -(grad_diff*grad_diff)/quad_coef;
00985 else
00986 obj_diff = -(grad_diff*grad_diff)/TAU;
00987
00988 if (obj_diff <= obj_diff_min)
00989 {
00990 Gmin_idx=j;
00991 obj_diff_min = obj_diff;
00992 }
00993 }
00994 }
00995 }
00996 }
00997
00998 gap=max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2);
00999 if(gap < eps)
01000 return 1;
01001
01002 if (y[Gmin_idx] == +1)
01003 out_i = Gmaxp_idx;
01004 else
01005 out_i = Gmaxn_idx;
01006 out_j = Gmin_idx;
01007
01008 return 0;
01009 }
01010
01011 bool Solver_NU::be_shrunk(
01012 int32_t i, float64_t Gmax1, float64_t Gmax2, float64_t Gmax3,
01013 float64_t Gmax4)
01014 {
01015 if(is_upper_bound(i))
01016 {
01017 if(y[i]==+1)
01018 return(-G[i] > Gmax1);
01019 else
01020 return(-G[i] > Gmax4);
01021 }
01022 else if(is_lower_bound(i))
01023 {
01024 if(y[i]==+1)
01025 return(G[i] > Gmax2);
01026 else
01027 return(G[i] > Gmax3);
01028 }
01029 else
01030 return(false);
01031 }
01032
01033 void Solver_NU::do_shrinking()
01034 {
01035 float64_t Gmax1 = -INF;
01036 float64_t Gmax2 = -INF;
01037 float64_t Gmax3 = -INF;
01038 float64_t Gmax4 = -INF;
01039
01040
01041 int32_t i;
01042 for(i=0;i<active_size;i++)
01043 {
01044 if(!is_upper_bound(i))
01045 {
01046 if(y[i]==+1)
01047 {
01048 if(-G[i] > Gmax1) Gmax1 = -G[i];
01049 }
01050 else if(-G[i] > Gmax4) Gmax4 = -G[i];
01051 }
01052 if(!is_lower_bound(i))
01053 {
01054 if(y[i]==+1)
01055 {
01056 if(G[i] > Gmax2) Gmax2 = G[i];
01057 }
01058 else if(G[i] > Gmax3) Gmax3 = G[i];
01059 }
01060 }
01061
01062 if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10)
01063 {
01064 unshrink = true;
01065 reconstruct_gradient();
01066 active_size = l;
01067 }
01068
01069 for(i=0;i<active_size;i++)
01070 if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
01071 {
01072 active_size--;
01073 while (active_size > i)
01074 {
01075 if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
01076 {
01077 swap_index(i,active_size);
01078 break;
01079 }
01080 active_size--;
01081 }
01082 }
01083 }
01084
01085 float64_t Solver_NU::calculate_rho()
01086 {
01087 int32_t nr_free1 = 0,nr_free2 = 0;
01088 float64_t ub1 = INF, ub2 = INF;
01089 float64_t lb1 = -INF, lb2 = -INF;
01090 float64_t sum_free1 = 0, sum_free2 = 0;
01091
01092 for(int32_t i=0; i<active_size; i++)
01093 {
01094 if(y[i]==+1)
01095 {
01096 if(is_upper_bound(i))
01097 lb1 = max(lb1,G[i]);
01098 else if(is_lower_bound(i))
01099 ub1 = min(ub1,G[i]);
01100 else
01101 {
01102 ++nr_free1;
01103 sum_free1 += G[i];
01104 }
01105 }
01106 else
01107 {
01108 if(is_upper_bound(i))
01109 lb2 = max(lb2,G[i]);
01110 else if(is_lower_bound(i))
01111 ub2 = min(ub2,G[i]);
01112 else
01113 {
01114 ++nr_free2;
01115 sum_free2 += G[i];
01116 }
01117 }
01118 }
01119
01120 float64_t r1,r2;
01121 if(nr_free1 > 0)
01122 r1 = sum_free1/nr_free1;
01123 else
01124 r1 = (ub1+lb1)/2;
01125
01126 if(nr_free2 > 0)
01127 r2 = sum_free2/nr_free2;
01128 else
01129 r2 = (ub2+lb2)/2;
01130
01131 si->r = (r1+r2)/2;
01132 return (r1-r2)/2;
01133 }
01134
01135
01136
01137
01138 class SVC_Q: public LibSVMKernel
01139 {
01140 public:
01141 SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
01142 :LibSVMKernel(prob.l, prob.x, param)
01143 {
01144 clone(y,y_,prob.l);
01145 cache = new Cache(prob.l,(int64_t)(param.cache_size*(1l<<20)));
01146 QD = new Qfloat[prob.l];
01147 for(int32_t i=0;i<prob.l;i++)
01148 QD[i]= (Qfloat)kernel_function(i,i);
01149 }
01150
01151 Qfloat *get_Q(int32_t i, int32_t len) const
01152 {
01153 Qfloat *data;
01154 int32_t start;
01155 if((start = cache->get_data(i,&data,len)) < len)
01156 {
01157 for(int32_t j=start;j<len;j++)
01158 data[j] = (Qfloat) y[i]*y[j]*kernel_function(i,j);
01159 }
01160 return data;
01161 }
01162
01163 Qfloat *get_QD() const
01164 {
01165 return QD;
01166 }
01167
01168 void swap_index(int32_t i, int32_t j) const
01169 {
01170 cache->swap_index(i,j);
01171 LibSVMKernel::swap_index(i,j);
01172 swap(y[i],y[j]);
01173 swap(QD[i],QD[j]);
01174 }
01175
01176 ~SVC_Q()
01177 {
01178 delete[] y;
01179 delete cache;
01180 delete[] QD;
01181 }
01182 private:
01183 schar *y;
01184 Cache *cache;
01185 Qfloat *QD;
01186 };
01187
01188 class ONE_CLASS_Q: public LibSVMKernel
01189 {
01190 public:
01191 ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
01192 :LibSVMKernel(prob.l, prob.x, param)
01193 {
01194 cache = new Cache(prob.l,(int64_t)(param.cache_size*(1l<<20)));
01195 QD = new Qfloat[prob.l];
01196 for(int32_t i=0;i<prob.l;i++)
01197 QD[i]= (Qfloat)kernel_function(i,i);
01198 }
01199
01200 Qfloat *get_Q(int32_t i, int32_t len) const
01201 {
01202 Qfloat *data;
01203 int32_t start;
01204 if((start = cache->get_data(i,&data,len)) < len)
01205 {
01206 for(int32_t j=start;j<len;j++)
01207 data[j] = (Qfloat) kernel_function(i,j);
01208 }
01209 return data;
01210 }
01211
01212 Qfloat *get_QD() const
01213 {
01214 return QD;
01215 }
01216
01217 void swap_index(int32_t i, int32_t j) const
01218 {
01219 cache->swap_index(i,j);
01220 LibSVMKernel::swap_index(i,j);
01221 swap(QD[i],QD[j]);
01222 }
01223
01224 ~ONE_CLASS_Q()
01225 {
01226 delete cache;
01227 delete[] QD;
01228 }
01229 private:
01230 Cache *cache;
01231 Qfloat *QD;
01232 };
01233
01234 class SVR_Q: public LibSVMKernel
01235 {
01236 public:
01237 SVR_Q(const svm_problem& prob, const svm_parameter& param)
01238 :LibSVMKernel(prob.l, prob.x, param)
01239 {
01240 l = prob.l;
01241 cache = new Cache(l,(int64_t)(param.cache_size*(1l<<20)));
01242 QD = new Qfloat[2*l];
01243 sign = new schar[2*l];
01244 index = new int32_t[2*l];
01245 for(int32_t k=0;k<l;k++)
01246 {
01247 sign[k] = 1;
01248 sign[k+l] = -1;
01249 index[k] = k;
01250 index[k+l] = k;
01251 QD[k]= (Qfloat)kernel_function(k,k);
01252 QD[k+l]=QD[k];
01253 }
01254 buffer[0] = new Qfloat[2*l];
01255 buffer[1] = new Qfloat[2*l];
01256 next_buffer = 0;
01257 }
01258
01259 void swap_index(int32_t i, int32_t j) const
01260 {
01261 swap(sign[i],sign[j]);
01262 swap(index[i],index[j]);
01263 swap(QD[i],QD[j]);
01264 }
01265
01266 Qfloat *get_Q(int32_t i, int32_t len) const
01267 {
01268 Qfloat *data;
01269 int32_t real_i = index[i];
01270 if(cache->get_data(real_i,&data,l) < l)
01271 {
01272 for(int32_t j=0;j<l;j++)
01273 data[j] = (Qfloat)kernel_function(real_i,j);
01274 }
01275
01276
01277 Qfloat *buf = buffer[next_buffer];
01278 next_buffer = 1 - next_buffer;
01279 schar si = sign[i];
01280 for(int32_t j=0;j<len;j++)
01281 buf[j] = si * sign[j] * data[index[j]];
01282 return buf;
01283 }
01284
01285 Qfloat *get_QD() const
01286 {
01287 return QD;
01288 }
01289
01290 ~SVR_Q()
01291 {
01292 delete cache;
01293 delete[] sign;
01294 delete[] index;
01295 delete[] buffer[0];
01296 delete[] buffer[1];
01297 delete[] QD;
01298 }
01299
01300 private:
01301 int32_t l;
01302 Cache *cache;
01303 schar *sign;
01304 int32_t *index;
01305 mutable int32_t next_buffer;
01306 Qfloat *buffer[2];
01307 Qfloat *QD;
01308 };
01309
01310
01311
01312
01313 static void solve_c_svc(
01314 const svm_problem *prob, const svm_parameter* param,
01315 float64_t *alpha, Solver::SolutionInfo* si, float64_t Cp, float64_t Cn)
01316 {
01317 int32_t l = prob->l;
01318 float64_t *minus_ones = new float64_t[l];
01319 schar *y = new schar[l];
01320
01321 int32_t i;
01322
01323 for(i=0;i<l;i++)
01324 {
01325 alpha[i] = 0;
01326 minus_ones[i] = -1;
01327 if(prob->y[i] > 0) y[i] = +1; else y[i]=-1;
01328 }
01329
01330 Solver s;
01331 s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
01332 alpha, Cp, Cn, param->eps, si, param->shrinking);
01333
01334 float64_t sum_alpha=0;
01335 for(i=0;i<l;i++)
01336 sum_alpha += alpha[i];
01337
01338 if (Cp==Cn)
01339 SG_SINFO("nu = %f\n", sum_alpha/(param->C*prob->l));
01340
01341 for(i=0;i<l;i++)
01342 alpha[i] *= y[i];
01343
01344 delete[] minus_ones;
01345 delete[] y;
01346 }
01347
01348 static void solve_nu_svc(
01349 const svm_problem *prob, const svm_parameter *param,
01350 float64_t *alpha, Solver::SolutionInfo* si)
01351 {
01352 int32_t i;
01353 int32_t l = prob->l;
01354 float64_t nu = param->nu;
01355
01356 schar *y = new schar[l];
01357
01358 for(i=0;i<l;i++)
01359 if(prob->y[i]>0)
01360 y[i] = +1;
01361 else
01362 y[i] = -1;
01363
01364 float64_t sum_pos = nu*l/2;
01365 float64_t sum_neg = nu*l/2;
01366
01367 for(i=0;i<l;i++)
01368 if(y[i] == +1)
01369 {
01370 alpha[i] = min(1.0,sum_pos);
01371 sum_pos -= alpha[i];
01372 }
01373 else
01374 {
01375 alpha[i] = min(1.0,sum_neg);
01376 sum_neg -= alpha[i];
01377 }
01378
01379 float64_t *zeros = new float64_t[l];
01380
01381 for(i=0;i<l;i++)
01382 zeros[i] = 0;
01383
01384 Solver_NU s;
01385 s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
01386 alpha, 1.0, 1.0, param->eps, si, param->shrinking);
01387 float64_t r = si->r;
01388
01389 SG_SINFO("C = %f\n",1/r);
01390
01391 for(i=0;i<l;i++)
01392 alpha[i] *= y[i]/r;
01393
01394 si->rho /= r;
01395 si->obj /= (r*r);
01396 si->upper_bound_p = 1/r;
01397 si->upper_bound_n = 1/r;
01398
01399 delete[] y;
01400 delete[] zeros;
01401 }
01402
01403 static void solve_one_class(
01404 const svm_problem *prob, const svm_parameter *param,
01405 float64_t *alpha, Solver::SolutionInfo* si)
01406 {
01407 int32_t l = prob->l;
01408 float64_t *zeros = new float64_t[l];
01409 schar *ones = new schar[l];
01410 int32_t i;
01411
01412 int32_t n = (int32_t)(param->nu*prob->l);
01413
01414 for(i=0;i<n;i++)
01415 alpha[i] = 1;
01416 if(n<prob->l)
01417 alpha[n] = param->nu * prob->l - n;
01418 for(i=n+1;i<l;i++)
01419 alpha[i] = 0;
01420
01421 for(i=0;i<l;i++)
01422 {
01423 zeros[i] = 0;
01424 ones[i] = 1;
01425 }
01426
01427 Solver s;
01428 s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
01429 alpha, 1.0, 1.0, param->eps, si, param->shrinking);
01430
01431 delete[] zeros;
01432 delete[] ones;
01433 }
01434
01435 static void solve_epsilon_svr(
01436 const svm_problem *prob, const svm_parameter *param,
01437 float64_t *alpha, Solver::SolutionInfo* si)
01438 {
01439 int32_t l = prob->l;
01440 float64_t *alpha2 = new float64_t[2*l];
01441 float64_t *linear_term = new float64_t[2*l];
01442 schar *y = new schar[2*l];
01443 int32_t i;
01444
01445 for(i=0;i<l;i++)
01446 {
01447 alpha2[i] = 0;
01448 linear_term[i] = param->p - prob->y[i];
01449 y[i] = 1;
01450
01451 alpha2[i+l] = 0;
01452 linear_term[i+l] = param->p + prob->y[i];
01453 y[i+l] = -1;
01454 }
01455
01456 Solver s;
01457 s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
01458 alpha2, param->C, param->C, param->eps, si, param->shrinking);
01459
01460 float64_t sum_alpha = 0;
01461 for(i=0;i<l;i++)
01462 {
01463 alpha[i] = alpha2[i] - alpha2[i+l];
01464 sum_alpha += fabs(alpha[i]);
01465 }
01466 SG_SINFO("nu = %f\n",sum_alpha/(param->C*l));
01467
01468 delete[] alpha2;
01469 delete[] linear_term;
01470 delete[] y;
01471 }
01472
01473 static void solve_nu_svr(
01474 const svm_problem *prob, const svm_parameter *param,
01475 float64_t *alpha, Solver::SolutionInfo* si)
01476 {
01477 int32_t l = prob->l;
01478 float64_t C = param->C;
01479 float64_t *alpha2 = new float64_t[2*l];
01480 float64_t *linear_term = new float64_t[2*l];
01481 schar *y = new schar[2*l];
01482 int32_t i;
01483
01484 float64_t sum = C * param->nu * l / 2;
01485 for(i=0;i<l;i++)
01486 {
01487 alpha2[i] = alpha2[i+l] = min(sum,C);
01488 sum -= alpha2[i];
01489
01490 linear_term[i] = - prob->y[i];
01491 y[i] = 1;
01492
01493 linear_term[i+l] = prob->y[i];
01494 y[i+l] = -1;
01495 }
01496
01497 Solver_NU s;
01498 s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
01499 alpha2, C, C, param->eps, si, param->shrinking);
01500
01501 SG_SINFO("epsilon = %f\n",-si->r);
01502
01503 for(i=0;i<l;i++)
01504 alpha[i] = alpha2[i] - alpha2[i+l];
01505
01506 delete[] alpha2;
01507 delete[] linear_term;
01508 delete[] y;
01509 }
01510
01511
01512
01513
01514 struct decision_function
01515 {
01516 float64_t *alpha;
01517 float64_t rho;
01518 float64_t objective;
01519 };
01520
01521 decision_function svm_train_one(
01522 const svm_problem *prob, const svm_parameter *param,
01523 float64_t Cp, float64_t Cn)
01524 {
01525 float64_t *alpha = Malloc(float64_t, prob->l);
01526 Solver::SolutionInfo si;
01527 switch(param->svm_type)
01528 {
01529 case C_SVC:
01530 solve_c_svc(prob,param,alpha,&si,Cp,Cn);
01531 break;
01532 case NU_SVC:
01533 solve_nu_svc(prob,param,alpha,&si);
01534 break;
01535 case ONE_CLASS:
01536 solve_one_class(prob,param,alpha,&si);
01537 break;
01538 case EPSILON_SVR:
01539 solve_epsilon_svr(prob,param,alpha,&si);
01540 break;
01541 case NU_SVR:
01542 solve_nu_svr(prob,param,alpha,&si);
01543 break;
01544 }
01545
01546 SG_SINFO("obj = %.16f, rho = %.16f\n",si.obj,si.rho);
01547
01548
01549 if (param->svm_type != ONE_CLASS)
01550 {
01551 int32_t nSV = 0;
01552 int32_t nBSV = 0;
01553 for(int32_t i=0;i<prob->l;i++)
01554 {
01555 if(fabs(alpha[i]) > 0)
01556 {
01557 ++nSV;
01558 if(prob->y[i] > 0)
01559 {
01560 if(fabs(alpha[i]) >= si.upper_bound_p)
01561 ++nBSV;
01562 }
01563 else
01564 {
01565 if(fabs(alpha[i]) >= si.upper_bound_n)
01566 ++nBSV;
01567 }
01568 }
01569 }
01570 SG_SINFO("nSV = %d, nBSV = %d\n",nSV,nBSV);
01571 }
01572
01573 decision_function f;
01574 f.alpha = alpha;
01575 f.rho = si.rho;
01576 f.objective=si.obj;
01577 return f;
01578 }
01579
01580
01581
01582
01583 void svm_group_classes(
01584 const svm_problem *prob, int32_t *nr_class_ret, int32_t **label_ret,
01585 int32_t **start_ret, int32_t **count_ret, int32_t *perm)
01586 {
01587 int32_t l = prob->l;
01588 int32_t max_nr_class = 16;
01589 int32_t nr_class = 0;
01590 int32_t *label = Malloc(int32_t, max_nr_class);
01591 int32_t *count = Malloc(int32_t, max_nr_class);
01592 int32_t *data_label = Malloc(int32_t, l);
01593 int32_t i;
01594
01595 for(i=0;i<l;i++)
01596 {
01597 int32_t this_label=(int32_t) prob->y[i];
01598 int32_t j;
01599 for(j=0;j<nr_class;j++)
01600 {
01601 if(this_label == label[j])
01602 {
01603 ++count[j];
01604 break;
01605 }
01606 }
01607 data_label[i] = j;
01608 if(j == nr_class)
01609 {
01610 if(nr_class == max_nr_class)
01611 {
01612 max_nr_class *= 2;
01613 label=(int32_t *) realloc(label,max_nr_class*sizeof(int32_t));
01614 count=(int32_t *) realloc(count,max_nr_class*sizeof(int32_t));
01615 }
01616 label[nr_class] = this_label;
01617 count[nr_class] = 1;
01618 ++nr_class;
01619 }
01620 }
01621
01622 int32_t *start = Malloc(int32_t, nr_class);
01623 start[0] = 0;
01624 for(i=1;i<nr_class;i++)
01625 start[i] = start[i-1]+count[i-1];
01626 for(i=0;i<l;i++)
01627 {
01628 perm[start[data_label[i]]] = i;
01629 ++start[data_label[i]];
01630 }
01631 start[0] = 0;
01632 for(i=1;i<nr_class;i++)
01633 start[i] = start[i-1]+count[i-1];
01634
01635 *nr_class_ret = nr_class;
01636 *label_ret = label;
01637 *start_ret = start;
01638 *count_ret = count;
01639 free(data_label);
01640 }
01641
01642
01643
01644
01645 svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
01646 {
01647 svm_model *model = Malloc(svm_model,1);
01648 model->param = *param;
01649 model->free_sv = 0;
01650
01651 if(param->svm_type == ONE_CLASS ||
01652 param->svm_type == EPSILON_SVR ||
01653 param->svm_type == NU_SVR)
01654 {
01655 SG_SINFO("training one class svm or doing epsilon sv regression\n");
01656
01657
01658 model->nr_class = 2;
01659 model->label = NULL;
01660 model->nSV = NULL;
01661 model->sv_coef = Malloc(float64_t *,1);
01662 decision_function f = svm_train_one(prob,param,0,0);
01663 model->rho = Malloc(float64_t, 1);
01664 model->rho[0] = f.rho;
01665 model->objective = f.objective;
01666
01667 int32_t nSV = 0;
01668 int32_t i;
01669 for(i=0;i<prob->l;i++)
01670 if(fabs(f.alpha[i]) > 0) ++nSV;
01671 model->l = nSV;
01672 model->SV = Malloc(svm_node *,nSV);
01673 model->sv_coef[0] = Malloc(float64_t, nSV);
01674 int32_t j = 0;
01675 for(i=0;i<prob->l;i++)
01676 if(fabs(f.alpha[i]) > 0)
01677 {
01678 model->SV[j] = prob->x[i];
01679 model->sv_coef[0][j] = f.alpha[i];
01680 ++j;
01681 }
01682
01683 free(f.alpha);
01684 }
01685 else
01686 {
01687
01688 int32_t l = prob->l;
01689 int32_t nr_class;
01690 int32_t *label = NULL;
01691 int32_t *start = NULL;
01692 int32_t *count = NULL;
01693 int32_t *perm = Malloc(int32_t, l);
01694
01695
01696 svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
01697 svm_node **x = Malloc(svm_node *,l);
01698 int32_t i;
01699 for(i=0;i<l;i++)
01700 x[i] = prob->x[perm[i]];
01701
01702
01703
01704 float64_t *weighted_C = Malloc(float64_t, nr_class);
01705 for(i=0;i<nr_class;i++)
01706 weighted_C[i] = param->C;
01707 for(i=0;i<param->nr_weight;i++)
01708 {
01709 int32_t j;
01710 for(j=0;j<nr_class;j++)
01711 if(param->weight_label[i] == label[j])
01712 break;
01713 if(j == nr_class)
01714 SG_SWARNING("warning: class label %d specified in weight is not found\n", param->weight_label[i]);
01715 else
01716 weighted_C[j] *= param->weight[i];
01717 }
01718
01719
01720
01721 bool *nonzero = Malloc(bool,l);
01722 for(i=0;i<l;i++)
01723 nonzero[i] = false;
01724 decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
01725
01726 int32_t p = 0;
01727 for(i=0;i<nr_class;i++)
01728 for(int32_t j=i+1;j<nr_class;j++)
01729 {
01730 svm_problem sub_prob;
01731 int32_t si = start[i], sj = start[j];
01732 int32_t ci = count[i], cj = count[j];
01733 sub_prob.l = ci+cj;
01734 sub_prob.x = Malloc(svm_node *,sub_prob.l);
01735 sub_prob.y = Malloc(float64_t, sub_prob.l+1);
01736
01737 int32_t k;
01738 for(k=0;k<ci;k++)
01739 {
01740 sub_prob.x[k] = x[si+k];
01741 sub_prob.y[k] = +1;
01742 }
01743 for(k=0;k<cj;k++)
01744 {
01745 sub_prob.x[ci+k] = x[sj+k];
01746 sub_prob.y[ci+k] = -1;
01747 }
01748 sub_prob.y[sub_prob.l]=-1;
01749
01750 f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
01751 for(k=0;k<ci;k++)
01752 if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
01753 nonzero[si+k] = true;
01754 for(k=0;k<cj;k++)
01755 if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
01756 nonzero[sj+k] = true;
01757 free(sub_prob.x);
01758 free(sub_prob.y);
01759 ++p;
01760 }
01761
01762
01763
01764 model->objective = f[0].objective;
01765 model->nr_class = nr_class;
01766
01767 model->label = Malloc(int32_t, nr_class);
01768 for(i=0;i<nr_class;i++)
01769 model->label[i] = label[i];
01770
01771 model->rho = Malloc(float64_t, nr_class*(nr_class-1)/2);
01772 for(i=0;i<nr_class*(nr_class-1)/2;i++)
01773 model->rho[i] = f[i].rho;
01774
01775 int32_t total_sv = 0;
01776 int32_t *nz_count = Malloc(int32_t, nr_class);
01777 model->nSV = Malloc(int32_t, nr_class);
01778 for(i=0;i<nr_class;i++)
01779 {
01780 int32_t nSV = 0;
01781 for(int32_t j=0;j<count[i];j++)
01782 if(nonzero[start[i]+j])
01783 {
01784 ++nSV;
01785 ++total_sv;
01786 }
01787 model->nSV[i] = nSV;
01788 nz_count[i] = nSV;
01789 }
01790
01791 SG_SINFO("Total nSV = %d\n",total_sv);
01792
01793 model->l = total_sv;
01794 model->SV = Malloc(svm_node *,total_sv);
01795 p = 0;
01796 for(i=0;i<l;i++)
01797 if(nonzero[i]) model->SV[p++] = x[i];
01798
01799 int32_t *nz_start = Malloc(int32_t, nr_class);
01800 nz_start[0] = 0;
01801 for(i=1;i<nr_class;i++)
01802 nz_start[i] = nz_start[i-1]+nz_count[i-1];
01803
01804 model->sv_coef = Malloc(float64_t *,nr_class-1);
01805 for(i=0;i<nr_class-1;i++)
01806 model->sv_coef[i] = Malloc(float64_t, total_sv);
01807
01808 p = 0;
01809 for(i=0;i<nr_class;i++)
01810 for(int32_t j=i+1;j<nr_class;j++)
01811 {
01812
01813
01814
01815
01816 int32_t si = start[i];
01817 int32_t sj = start[j];
01818 int32_t ci = count[i];
01819 int32_t cj = count[j];
01820
01821 int32_t q = nz_start[i];
01822 int32_t k;
01823 for(k=0;k<ci;k++)
01824 if(nonzero[si+k])
01825 model->sv_coef[j-1][q++] = f[p].alpha[k];
01826 q = nz_start[j];
01827 for(k=0;k<cj;k++)
01828 if(nonzero[sj+k])
01829 model->sv_coef[i][q++] = f[p].alpha[ci+k];
01830 ++p;
01831 }
01832
01833 free(label);
01834 free(count);
01835 free(perm);
01836 free(start);
01837 free(x);
01838 free(weighted_C);
01839 free(nonzero);
01840 for(i=0;i<nr_class*(nr_class-1)/2;i++)
01841 free(f[i].alpha);
01842 free(f);
01843 free(nz_count);
01844 free(nz_start);
01845 }
01846 return model;
01847 }
01848
01849 const char *svm_type_table[] =
01850 {
01851 "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
01852 };
01853
01854 const char *kernel_type_table[]=
01855 {
01856 "linear","polynomial","rbf","sigmoid","precomputed",NULL
01857 };
01858
01859 void svm_destroy_model(svm_model* model)
01860 {
01861 if(model->free_sv && model->l > 0)
01862 free((void *)(model->SV[0]));
01863 for(int32_t i=0;i<model->nr_class-1;i++)
01864 free(model->sv_coef[i]);
01865 free(model->SV);
01866 free(model->sv_coef);
01867 free(model->rho);
01868 free(model->label);
01869 free(model->nSV);
01870 free(model);
01871 }
01872
01873 void svm_destroy_param(svm_parameter* param)
01874 {
01875 free(param->weight_label);
01876 free(param->weight);
01877 }
01878
01879 const char *svm_check_parameter(
01880 const svm_problem *prob, const svm_parameter *param)
01881 {
01882
01883
01884 int32_t svm_type = param->svm_type;
01885 if(svm_type != C_SVC &&
01886 svm_type != NU_SVC &&
01887 svm_type != ONE_CLASS &&
01888 svm_type != EPSILON_SVR &&
01889 svm_type != NU_SVR)
01890 return "unknown svm type";
01891
01892
01893
01894 int32_t kernel_type = param->kernel_type;
01895 if(kernel_type != LINEAR &&
01896 kernel_type != POLY &&
01897 kernel_type != RBF &&
01898 kernel_type != SIGMOID &&
01899 kernel_type != PRECOMPUTED)
01900 return "unknown kernel type";
01901
01902 if(param->degree < 0)
01903 return "degree of polynomial kernel < 0";
01904
01905
01906
01907 if(param->cache_size <= 0)
01908 return "cache_size <= 0";
01909
01910 if(param->eps <= 0)
01911 return "eps <= 0";
01912
01913 if(svm_type == C_SVC ||
01914 svm_type == EPSILON_SVR ||
01915 svm_type == NU_SVR)
01916 if(param->C <= 0)
01917 return "C <= 0";
01918
01919 if(svm_type == NU_SVC ||
01920 svm_type == ONE_CLASS ||
01921 svm_type == NU_SVR)
01922 if(param->nu <= 0 || param->nu > 1)
01923 return "nu <= 0 or nu > 1";
01924
01925 if(svm_type == EPSILON_SVR)
01926 if(param->p < 0)
01927 return "p < 0";
01928
01929 if(param->shrinking != 0 &&
01930 param->shrinking != 1)
01931 return "shrinking != 0 and shrinking != 1";
01932
01933
01934
01935
01936 if(svm_type == NU_SVC)
01937 {
01938 int32_t l = prob->l;
01939 int32_t max_nr_class = 16;
01940 int32_t nr_class = 0;
01941 int32_t *label = Malloc(int32_t, max_nr_class);
01942 int32_t *count = Malloc(int32_t, max_nr_class);
01943
01944 int32_t i;
01945 for(i=0;i<l;i++)
01946 {
01947 int32_t this_label = (int32_t) prob->y[i];
01948 int32_t j;
01949 for(j=0;j<nr_class;j++)
01950 if(this_label == label[j])
01951 {
01952 ++count[j];
01953 break;
01954 }
01955 if(j == nr_class)
01956 {
01957 if(nr_class == max_nr_class)
01958 {
01959 max_nr_class *= 2;
01960 label=(int32_t *) realloc(label,
01961 max_nr_class*sizeof(int32_t));
01962 count=(int32_t *) realloc(count,
01963 max_nr_class*sizeof(int32_t));
01964 }
01965 label[nr_class] = this_label;
01966 count[nr_class] = 1;
01967 ++nr_class;
01968 }
01969 }
01970
01971 for(i=0;i<nr_class;i++)
01972 {
01973 int32_t n1 = count[i];
01974 for(int32_t j=i+1;j<nr_class;j++)
01975 {
01976 int32_t n2 = count[j];
01977 if(param->nu*(n1+n2)/2 > min(n1,n2))
01978 {
01979 free(label);
01980 free(count);
01981 return "specified nu is infeasible";
01982 }
01983 }
01984 }
01985 free(label);
01986 free(count);
01987 }
01988
01989 return NULL;
01990 }
01991 #endif // DOXYGEN_SHOULD_SKIP_THIS