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