00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <ctype.h>
00005 #include <float.h>
00006 #include <string.h>
00007 #include <stdarg.h>
00008 #include "svm.h"
00009 int libsvm_version = LIBSVM_VERSION;
00010 typedef float Qfloat;
00011 typedef signed char schar;
00012 #ifndef min
00013 template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
00014 #endif
00015 #ifndef max
00016 template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
00017 #endif
00018 template <class T> static inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
00019 template <class S, class T> static inline void clone(T*& dst, S* src, int n)
00020 {
00021 dst = new T[n];
00022 memcpy((void *)dst,(void *)src,sizeof(T)*n);
00023 }
00024 static inline double powi(double base, int times)
00025 {
00026 double tmp = base, ret = 1.0;
00027
00028 for(int t=times; t>0; t/=2)
00029 {
00030 if(t%2==1) ret*=tmp;
00031 tmp = tmp * tmp;
00032 }
00033 return ret;
00034 }
00035 #define INF HUGE_VAL
00036 #define TAU 1e-12
00037 #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
00038
00039 static void print_string_stdout(const char *s)
00040 {
00041 fputs(s,stdout);
00042 fflush(stdout);
00043 }
00044 static void (*svm_print_string) (const char *) = &print_string_stdout;
00045 #if 1
00046 static void info(const char *fmt,...)
00047 {
00048 char buf[BUFSIZ];
00049 va_list ap;
00050 va_start(ap,fmt);
00051 vsprintf(buf,fmt,ap);
00052 va_end(ap);
00053 (*svm_print_string)(buf);
00054 }
00055 #else
00056 static void info(const char *fmt,...) {}
00057 #endif
00058
00059
00060
00061
00062
00063
00064
00065 class Cache
00066 {
00067 public:
00068 Cache(int l,long int size);
00069 ~Cache();
00070
00071
00072
00073
00074 int get_data(const int index, Qfloat **data, int len);
00075 void swap_index(int i, int j);
00076 private:
00077 int l;
00078 long int size;
00079 struct head_t
00080 {
00081 head_t *prev, *next;
00082 Qfloat *data;
00083 int len;
00084 };
00085
00086 head_t *head;
00087 head_t lru_head;
00088 void lru_delete(head_t *h);
00089 void lru_insert(head_t *h);
00090 };
00091
00092 Cache::Cache(int l_,long int size_):l(l_),size(size_)
00093 {
00094 head = (head_t *)calloc(l,sizeof(head_t));
00095 size /= sizeof(Qfloat);
00096 size -= l * sizeof(head_t) / sizeof(Qfloat);
00097 size = max(size, 2 * (long int) l);
00098 lru_head.next = lru_head.prev = &lru_head;
00099 }
00100
00101 Cache::~Cache()
00102 {
00103 for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
00104 free(h->data);
00105 free(head);
00106 }
00107
00108 void Cache::lru_delete(head_t *h)
00109 {
00110
00111 h->prev->next = h->next;
00112 h->next->prev = h->prev;
00113 }
00114
00115 void Cache::lru_insert(head_t *h)
00116 {
00117
00118 h->next = &lru_head;
00119 h->prev = lru_head.prev;
00120 h->prev->next = h;
00121 h->next->prev = h;
00122 }
00123
00124 int Cache::get_data(const int index, Qfloat **data, int len)
00125 {
00126 head_t *h = &head[index];
00127 if(h->len) lru_delete(h);
00128 int more = len - h->len;
00129
00130 if(more > 0)
00131 {
00132
00133 while(size < more)
00134 {
00135 head_t *old = lru_head.next;
00136 lru_delete(old);
00137 free(old->data);
00138 size += old->len;
00139 old->data = 0;
00140 old->len = 0;
00141 }
00142
00143
00144 h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
00145 size -= more;
00146 swap(h->len,len);
00147 }
00148
00149 lru_insert(h);
00150 *data = h->data;
00151 return len;
00152 }
00153
00154 void Cache::swap_index(int i, int j)
00155 {
00156 if(i==j) return;
00157
00158 if(head[i].len) lru_delete(&head[i]);
00159 if(head[j].len) lru_delete(&head[j]);
00160 swap(head[i].data,head[j].data);
00161 swap(head[i].len,head[j].len);
00162 if(head[i].len) lru_insert(&head[i]);
00163 if(head[j].len) lru_insert(&head[j]);
00164
00165 if(i>j) swap(i,j);
00166 for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
00167 {
00168 if(h->len > i)
00169 {
00170 if(h->len > j)
00171 swap(h->data[i],h->data[j]);
00172 else
00173 {
00174
00175 lru_delete(h);
00176 free(h->data);
00177 size += h->len;
00178 h->data = 0;
00179 h->len = 0;
00180 }
00181 }
00182 }
00183 }
00184
00185
00186
00187
00188
00189
00190
00191
00192 class QMatrix {
00193 public:
00194 virtual Qfloat *get_Q(int column, int len) const = 0;
00195 virtual double *get_QD() const = 0;
00196 virtual void swap_index(int i, int j) const = 0;
00197 virtual ~QMatrix() {}
00198 };
00199
00200 class Kernel: public QMatrix {
00201 public:
00202 Kernel(int l, svm_node * const * x, const svm_parameter& param);
00203 virtual ~Kernel();
00204
00205 static double k_function(const svm_node *x, const svm_node *y,
00206 const svm_parameter& param);
00207 virtual Qfloat *get_Q(int column, int len) const = 0;
00208 virtual double *get_QD() const = 0;
00209 virtual void swap_index(int i, int j) const
00210 {
00211 swap(x[i],x[j]);
00212 if(x_square) swap(x_square[i],x_square[j]);
00213 }
00214 protected:
00215
00216 double (Kernel::*kernel_function)(int i, int j) const;
00217
00218 private:
00219 const svm_node **x;
00220 double *x_square;
00221
00222
00223 const int kernel_type;
00224 const int degree;
00225 const double gamma;
00226 const double coef0;
00227
00228 static double dot(const svm_node *px, const svm_node *py);
00229 double kernel_linear(int i, int j) const
00230 {
00231 return dot(x[i],x[j]);
00232 }
00233 double kernel_poly(int i, int j) const
00234 {
00235 return powi(gamma*dot(x[i],x[j])+coef0,degree);
00236 }
00237 double kernel_rbf(int i, int j) const
00238 {
00239 return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
00240 }
00241 double kernel_sigmoid(int i, int j) const
00242 {
00243 return tanh(gamma*dot(x[i],x[j])+coef0);
00244 }
00245 double kernel_precomputed(int i, int j) const
00246 {
00247 return x[i][(int)(x[j][0].value)].value;
00248 }
00249 };
00250
00251 Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
00252 :kernel_type(param.kernel_type), degree(param.degree),
00253 gamma(param.gamma), coef0(param.coef0)
00254 {
00255 switch(kernel_type)
00256 {
00257 case LINEAR:
00258 kernel_function = &Kernel::kernel_linear;
00259 break;
00260 case POLY:
00261 kernel_function = &Kernel::kernel_poly;
00262 break;
00263 case RBF:
00264 kernel_function = &Kernel::kernel_rbf;
00265 break;
00266 case SIGMOID:
00267 kernel_function = &Kernel::kernel_sigmoid;
00268 break;
00269 case PRECOMPUTED:
00270 kernel_function = &Kernel::kernel_precomputed;
00271 break;
00272 }
00273
00274 clone(x,x_,l);
00275
00276 if(kernel_type == RBF)
00277 {
00278 x_square = new double[l];
00279 for(int i=0;i<l;i++)
00280 x_square[i] = dot(x[i],x[i]);
00281 }
00282 else
00283 x_square = 0;
00284 }
00285
00286 Kernel::~Kernel()
00287 {
00288 delete[] x;
00289 delete[] x_square;
00290 }
00291
00292 double Kernel::dot(const svm_node *px, const svm_node *py)
00293 {
00294 double sum = 0;
00295 while(px->index != -1 && py->index != -1)
00296 {
00297 if(px->index == py->index)
00298 {
00299 sum += px->value * py->value;
00300 ++px;
00301 ++py;
00302 }
00303 else
00304 {
00305 if(px->index > py->index)
00306 ++py;
00307 else
00308 ++px;
00309 }
00310 }
00311 return sum;
00312 }
00313
00314 double Kernel::k_function(const svm_node *x, const svm_node *y,
00315 const svm_parameter& param)
00316 {
00317 switch(param.kernel_type)
00318 {
00319 case LINEAR:
00320 return dot(x,y);
00321 case POLY:
00322 return powi(param.gamma*dot(x,y)+param.coef0,param.degree);
00323 case RBF:
00324 {
00325 double sum = 0;
00326 while(x->index != -1 && y->index !=-1)
00327 {
00328 if(x->index == y->index)
00329 {
00330 double d = x->value - y->value;
00331 sum += d*d;
00332 ++x;
00333 ++y;
00334 }
00335 else
00336 {
00337 if(x->index > y->index)
00338 {
00339 sum += y->value * y->value;
00340 ++y;
00341 }
00342 else
00343 {
00344 sum += x->value * x->value;
00345 ++x;
00346 }
00347 }
00348 }
00349
00350 while(x->index != -1)
00351 {
00352 sum += x->value * x->value;
00353 ++x;
00354 }
00355
00356 while(y->index != -1)
00357 {
00358 sum += y->value * y->value;
00359 ++y;
00360 }
00361
00362 return exp(-param.gamma*sum);
00363 }
00364 case SIGMOID:
00365 return tanh(param.gamma*dot(x,y)+param.coef0);
00366 case PRECOMPUTED:
00367 return x[(int)(y->value)].value;
00368 default:
00369 return 0;
00370 }
00371 }
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 class Solver {
00392 public:
00393 Solver() {};
00394 virtual ~Solver() {};
00395
00396 struct SolutionInfo {
00397 double obj;
00398 double rho;
00399 double upper_bound_p;
00400 double upper_bound_n;
00401 double r;
00402 };
00403
00404 void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
00405 double *alpha_, double Cp, double Cn, double eps,
00406 SolutionInfo* si, int shrinking);
00407 protected:
00408 int active_size;
00409 schar *y;
00410 double *G;
00411 enum { LOWER_BOUND, UPPER_BOUND, FREE };
00412 char *alpha_status;
00413 double *alpha;
00414 const QMatrix *Q;
00415 const double *QD;
00416 double eps;
00417 double Cp,Cn;
00418 double *p;
00419 int *active_set;
00420 double *G_bar;
00421 int l;
00422 bool unshrink;
00423
00424 double get_C(int i)
00425 {
00426 return (y[i] > 0)? Cp : Cn;
00427 }
00428 void update_alpha_status(int i)
00429 {
00430 if(alpha[i] >= get_C(i))
00431 alpha_status[i] = UPPER_BOUND;
00432 else if(alpha[i] <= 0)
00433 alpha_status[i] = LOWER_BOUND;
00434 else alpha_status[i] = FREE;
00435 }
00436 bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
00437 bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
00438 bool is_free(int i) { return alpha_status[i] == FREE; }
00439 void swap_index(int i, int j);
00440 void reconstruct_gradient();
00441 virtual int select_working_set(int &i, int &j);
00442 virtual double calculate_rho();
00443 virtual void do_shrinking();
00444 private:
00445 bool be_shrunk(int i, double Gmax1, double Gmax2);
00446 };
00447
00448 void Solver::swap_index(int i, int j)
00449 {
00450 Q->swap_index(i,j);
00451 swap(y[i],y[j]);
00452 swap(G[i],G[j]);
00453 swap(alpha_status[i],alpha_status[j]);
00454 swap(alpha[i],alpha[j]);
00455 swap(p[i],p[j]);
00456 swap(active_set[i],active_set[j]);
00457 swap(G_bar[i],G_bar[j]);
00458 }
00459
00460 void Solver::reconstruct_gradient()
00461 {
00462
00463
00464 if(active_size == l) return;
00465
00466 int i,j;
00467 int nr_free = 0;
00468
00469 for(j=active_size;j<l;j++)
00470 G[j] = G_bar[j] + p[j];
00471
00472 for(j=0;j<active_size;j++)
00473 if(is_free(j))
00474 nr_free++;
00475
00476 if(2*nr_free < active_size)
00477 info("\nWarning: using -h 0 may be faster\n");
00478
00479 if (nr_free*l > 2*active_size*(l-active_size))
00480 {
00481 for(i=active_size;i<l;i++)
00482 {
00483 const Qfloat *Q_i = Q->get_Q(i,active_size);
00484 for(j=0;j<active_size;j++)
00485 if(is_free(j))
00486 G[i] += alpha[j] * Q_i[j];
00487 }
00488 }
00489 else
00490 {
00491 for(i=0;i<active_size;i++)
00492 if(is_free(i))
00493 {
00494 const Qfloat *Q_i = Q->get_Q(i,l);
00495 double alpha_i = alpha[i];
00496 for(j=active_size;j<l;j++)
00497 G[j] += alpha_i * Q_i[j];
00498 }
00499 }
00500 }
00501
00502 void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
00503 double *alpha_, double Cp, double Cn, double eps,
00504 SolutionInfo* si, int shrinking)
00505 {
00506 this->l = l;
00507 this->Q = &Q;
00508 QD=Q.get_QD();
00509 clone(p, p_,l);
00510 clone(y, y_,l);
00511 clone(alpha,alpha_,l);
00512 this->Cp = Cp;
00513 this->Cn = Cn;
00514 this->eps = eps;
00515 unshrink = false;
00516
00517
00518 {
00519 alpha_status = new char[l];
00520 for(int i=0;i<l;i++)
00521 update_alpha_status(i);
00522 }
00523
00524
00525 {
00526 active_set = new int[l];
00527 for(int i=0;i<l;i++)
00528 active_set[i] = i;
00529 active_size = l;
00530 }
00531
00532
00533 {
00534 G = new double[l];
00535 G_bar = new double[l];
00536 int i;
00537 for(i=0;i<l;i++)
00538 {
00539 G[i] = p[i];
00540 G_bar[i] = 0;
00541 }
00542 for(i=0;i<l;i++)
00543 if(!is_lower_bound(i))
00544 {
00545 const Qfloat *Q_i = Q.get_Q(i,l);
00546 double alpha_i = alpha[i];
00547 int j;
00548 for(j=0;j<l;j++)
00549 G[j] += alpha_i*Q_i[j];
00550 if(is_upper_bound(i))
00551 for(j=0;j<l;j++)
00552 G_bar[j] += get_C(i) * Q_i[j];
00553 }
00554 }
00555
00556
00557
00558 int iter = 0;
00559 int counter = min(l,1000)+1;
00560
00561 while(1)
00562 {
00563
00564
00565 if(--counter == 0)
00566 {
00567 counter = min(l,1000);
00568 if(shrinking) do_shrinking();
00569 info(".");
00570 }
00571
00572 int i,j;
00573 if(select_working_set(i,j)!=0)
00574 {
00575
00576 reconstruct_gradient();
00577
00578 active_size = l;
00579 info("*");
00580 if(select_working_set(i,j)!=0)
00581 break;
00582 else
00583 counter = 1;
00584 }
00585
00586 ++iter;
00587
00588
00589
00590 const Qfloat *Q_i = Q.get_Q(i,active_size);
00591 const Qfloat *Q_j = Q.get_Q(j,active_size);
00592
00593 double C_i = get_C(i);
00594 double C_j = get_C(j);
00595
00596 double old_alpha_i = alpha[i];
00597 double old_alpha_j = alpha[j];
00598
00599 if(y[i]!=y[j])
00600 {
00601 double quad_coef = QD[i]+QD[j]+2*Q_i[j];
00602 if (quad_coef <= 0)
00603 quad_coef = TAU;
00604 double delta = (-G[i]-G[j])/quad_coef;
00605 double diff = alpha[i] - alpha[j];
00606 alpha[i] += delta;
00607 alpha[j] += delta;
00608
00609 if(diff > 0)
00610 {
00611 if(alpha[j] < 0)
00612 {
00613 alpha[j] = 0;
00614 alpha[i] = diff;
00615 }
00616 }
00617 else
00618 {
00619 if(alpha[i] < 0)
00620 {
00621 alpha[i] = 0;
00622 alpha[j] = -diff;
00623 }
00624 }
00625 if(diff > C_i - C_j)
00626 {
00627 if(alpha[i] > C_i)
00628 {
00629 alpha[i] = C_i;
00630 alpha[j] = C_i - diff;
00631 }
00632 }
00633 else
00634 {
00635 if(alpha[j] > C_j)
00636 {
00637 alpha[j] = C_j;
00638 alpha[i] = C_j + diff;
00639 }
00640 }
00641 }
00642 else
00643 {
00644 double quad_coef = QD[i]+QD[j]-2*Q_i[j];
00645 if (quad_coef <= 0)
00646 quad_coef = TAU;
00647 double delta = (G[i]-G[j])/quad_coef;
00648 double sum = alpha[i] + alpha[j];
00649 alpha[i] -= delta;
00650 alpha[j] += delta;
00651
00652 if(sum > C_i)
00653 {
00654 if(alpha[i] > C_i)
00655 {
00656 alpha[i] = C_i;
00657 alpha[j] = sum - C_i;
00658 }
00659 }
00660 else
00661 {
00662 if(alpha[j] < 0)
00663 {
00664 alpha[j] = 0;
00665 alpha[i] = sum;
00666 }
00667 }
00668 if(sum > C_j)
00669 {
00670 if(alpha[j] > C_j)
00671 {
00672 alpha[j] = C_j;
00673 alpha[i] = sum - C_j;
00674 }
00675 }
00676 else
00677 {
00678 if(alpha[i] < 0)
00679 {
00680 alpha[i] = 0;
00681 alpha[j] = sum;
00682 }
00683 }
00684 }
00685
00686
00687
00688 double delta_alpha_i = alpha[i] - old_alpha_i;
00689 double delta_alpha_j = alpha[j] - old_alpha_j;
00690
00691 for(int k=0;k<active_size;k++)
00692 {
00693 G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
00694 }
00695
00696
00697
00698 {
00699 bool ui = is_upper_bound(i);
00700 bool uj = is_upper_bound(j);
00701 update_alpha_status(i);
00702 update_alpha_status(j);
00703 int k;
00704 if(ui != is_upper_bound(i))
00705 {
00706 Q_i = Q.get_Q(i,l);
00707 if(ui)
00708 for(k=0;k<l;k++)
00709 G_bar[k] -= C_i * Q_i[k];
00710 else
00711 for(k=0;k<l;k++)
00712 G_bar[k] += C_i * Q_i[k];
00713 }
00714
00715 if(uj != is_upper_bound(j))
00716 {
00717 Q_j = Q.get_Q(j,l);
00718 if(uj)
00719 for(k=0;k<l;k++)
00720 G_bar[k] -= C_j * Q_j[k];
00721 else
00722 for(k=0;k<l;k++)
00723 G_bar[k] += C_j * Q_j[k];
00724 }
00725 }
00726 }
00727
00728
00729
00730 si->rho = calculate_rho();
00731
00732
00733 {
00734 double v = 0;
00735 int i;
00736 for(i=0;i<l;i++)
00737 v += alpha[i] * (G[i] + p[i]);
00738
00739 si->obj = v/2;
00740 }
00741
00742
00743 {
00744 for(int i=0;i<l;i++)
00745 alpha_[active_set[i]] = alpha[i];
00746 }
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756 si->upper_bound_p = Cp;
00757 si->upper_bound_n = Cn;
00758
00759 info("\noptimization finished, #iter = %d\n",iter);
00760
00761 delete[] p;
00762 delete[] y;
00763 delete[] alpha;
00764 delete[] alpha_status;
00765 delete[] active_set;
00766 delete[] G;
00767 delete[] G_bar;
00768 }
00769
00770
00771 int Solver::select_working_set(int &out_i, int &out_j)
00772 {
00773
00774
00775
00776
00777
00778
00779 double Gmax = -INF;
00780 double Gmax2 = -INF;
00781 int Gmax_idx = -1;
00782 int Gmin_idx = -1;
00783 double obj_diff_min = INF;
00784
00785 for(int t=0;t<active_size;t++)
00786 if(y[t]==+1)
00787 {
00788 if(!is_upper_bound(t))
00789 if(-G[t] >= Gmax)
00790 {
00791 Gmax = -G[t];
00792 Gmax_idx = t;
00793 }
00794 }
00795 else
00796 {
00797 if(!is_lower_bound(t))
00798 if(G[t] >= Gmax)
00799 {
00800 Gmax = G[t];
00801 Gmax_idx = t;
00802 }
00803 }
00804
00805 int i = Gmax_idx;
00806 const Qfloat *Q_i = NULL;
00807 if(i != -1)
00808 Q_i = Q->get_Q(i,active_size);
00809
00810 for(int j=0;j<active_size;j++)
00811 {
00812 if(y[j]==+1)
00813 {
00814 if (!is_lower_bound(j))
00815 {
00816 double grad_diff=Gmax+G[j];
00817 if (G[j] >= Gmax2)
00818 Gmax2 = G[j];
00819 if (grad_diff > 0)
00820 {
00821 double obj_diff;
00822 double quad_coef = QD[i]+QD[j]-2.0*y[i]*Q_i[j];
00823 if (quad_coef > 0)
00824 obj_diff = -(grad_diff*grad_diff)/quad_coef;
00825 else
00826 obj_diff = -(grad_diff*grad_diff)/TAU;
00827
00828 if (obj_diff <= obj_diff_min)
00829 {
00830 Gmin_idx=j;
00831 obj_diff_min = obj_diff;
00832 }
00833 }
00834 }
00835 }
00836 else
00837 {
00838 if (!is_upper_bound(j))
00839 {
00840 double grad_diff= Gmax-G[j];
00841 if (-G[j] >= Gmax2)
00842 Gmax2 = -G[j];
00843 if (grad_diff > 0)
00844 {
00845 double obj_diff;
00846 double quad_coef = QD[i]+QD[j]+2.0*y[i]*Q_i[j];
00847 if (quad_coef > 0)
00848 obj_diff = -(grad_diff*grad_diff)/quad_coef;
00849 else
00850 obj_diff = -(grad_diff*grad_diff)/TAU;
00851
00852 if (obj_diff <= obj_diff_min)
00853 {
00854 Gmin_idx=j;
00855 obj_diff_min = obj_diff;
00856 }
00857 }
00858 }
00859 }
00860 }
00861
00862 if(Gmax+Gmax2 < eps)
00863 return 1;
00864
00865 out_i = Gmax_idx;
00866 out_j = Gmin_idx;
00867 return 0;
00868 }
00869
00870 bool Solver::be_shrunk(int i, double Gmax1, double Gmax2)
00871 {
00872 if(is_upper_bound(i))
00873 {
00874 if(y[i]==+1)
00875 return(-G[i] > Gmax1);
00876 else
00877 return(-G[i] > Gmax2);
00878 }
00879 else if(is_lower_bound(i))
00880 {
00881 if(y[i]==+1)
00882 return(G[i] > Gmax2);
00883 else
00884 return(G[i] > Gmax1);
00885 }
00886 else
00887 return(false);
00888 }
00889
00890 void Solver::do_shrinking()
00891 {
00892 int i;
00893 double Gmax1 = -INF;
00894 double Gmax2 = -INF;
00895
00896
00897 for(i=0;i<active_size;i++)
00898 {
00899 if(y[i]==+1)
00900 {
00901 if(!is_upper_bound(i))
00902 {
00903 if(-G[i] >= Gmax1)
00904 Gmax1 = -G[i];
00905 }
00906 if(!is_lower_bound(i))
00907 {
00908 if(G[i] >= Gmax2)
00909 Gmax2 = G[i];
00910 }
00911 }
00912 else
00913 {
00914 if(!is_upper_bound(i))
00915 {
00916 if(-G[i] >= Gmax2)
00917 Gmax2 = -G[i];
00918 }
00919 if(!is_lower_bound(i))
00920 {
00921 if(G[i] >= Gmax1)
00922 Gmax1 = G[i];
00923 }
00924 }
00925 }
00926
00927 if(unshrink == false && Gmax1 + Gmax2 <= eps*10)
00928 {
00929 unshrink = true;
00930 reconstruct_gradient();
00931 active_size = l;
00932 info("*");
00933 }
00934
00935 for(i=0;i<active_size;i++)
00936 if (be_shrunk(i, Gmax1, Gmax2))
00937 {
00938 active_size--;
00939 while (active_size > i)
00940 {
00941 if (!be_shrunk(active_size, Gmax1, Gmax2))
00942 {
00943 swap_index(i,active_size);
00944 break;
00945 }
00946 active_size--;
00947 }
00948 }
00949 }
00950
00951 double Solver::calculate_rho()
00952 {
00953 double r;
00954 int nr_free = 0;
00955 double ub = INF, lb = -INF, sum_free = 0;
00956 for(int i=0;i<active_size;i++)
00957 {
00958 double yG = y[i]*G[i];
00959
00960 if(is_upper_bound(i))
00961 {
00962 if(y[i]==-1)
00963 ub = min(ub,yG);
00964 else
00965 lb = max(lb,yG);
00966 }
00967 else if(is_lower_bound(i))
00968 {
00969 if(y[i]==+1)
00970 ub = min(ub,yG);
00971 else
00972 lb = max(lb,yG);
00973 }
00974 else
00975 {
00976 ++nr_free;
00977 sum_free += yG;
00978 }
00979 }
00980
00981 if(nr_free>0)
00982 r = sum_free/nr_free;
00983 else
00984 r = (ub+lb)/2;
00985
00986 return r;
00987 }
00988
00989
00990
00991
00992
00993
00994 class Solver_NU : public Solver
00995 {
00996 public:
00997 Solver_NU() {}
00998 void Solve(int l, const QMatrix& Q, const double *p, const schar *y,
00999 double *alpha, double Cp, double Cn, double eps,
01000 SolutionInfo* si, int shrinking)
01001 {
01002 this->si = si;
01003 Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking);
01004 }
01005 private:
01006 SolutionInfo *si;
01007 int select_working_set(int &i, int &j);
01008 double calculate_rho();
01009 bool be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4);
01010 void do_shrinking();
01011 };
01012
01013
01014 int Solver_NU::select_working_set(int &out_i, int &out_j)
01015 {
01016
01017
01018
01019
01020
01021
01022 double Gmaxp = -INF;
01023 double Gmaxp2 = -INF;
01024 int Gmaxp_idx = -1;
01025
01026 double Gmaxn = -INF;
01027 double Gmaxn2 = -INF;
01028 int Gmaxn_idx = -1;
01029
01030 int Gmin_idx = -1;
01031 double obj_diff_min = INF;
01032
01033 for(int t=0;t<active_size;t++)
01034 if(y[t]==+1)
01035 {
01036 if(!is_upper_bound(t))
01037 if(-G[t] >= Gmaxp)
01038 {
01039 Gmaxp = -G[t];
01040 Gmaxp_idx = t;
01041 }
01042 }
01043 else
01044 {
01045 if(!is_lower_bound(t))
01046 if(G[t] >= Gmaxn)
01047 {
01048 Gmaxn = G[t];
01049 Gmaxn_idx = t;
01050 }
01051 }
01052
01053 int ip = Gmaxp_idx;
01054 int in = Gmaxn_idx;
01055 const Qfloat *Q_ip = NULL;
01056 const Qfloat *Q_in = NULL;
01057 if(ip != -1)
01058 Q_ip = Q->get_Q(ip,active_size);
01059 if(in != -1)
01060 Q_in = Q->get_Q(in,active_size);
01061
01062 for(int j=0;j<active_size;j++)
01063 {
01064 if(y[j]==+1)
01065 {
01066 if (!is_lower_bound(j))
01067 {
01068 double grad_diff=Gmaxp+G[j];
01069 if (G[j] >= Gmaxp2)
01070 Gmaxp2 = G[j];
01071 if (grad_diff > 0)
01072 {
01073 double obj_diff;
01074 double quad_coef = QD[ip]+QD[j]-2*Q_ip[j];
01075 if (quad_coef > 0)
01076 obj_diff = -(grad_diff*grad_diff)/quad_coef;
01077 else
01078 obj_diff = -(grad_diff*grad_diff)/TAU;
01079
01080 if (obj_diff <= obj_diff_min)
01081 {
01082 Gmin_idx=j;
01083 obj_diff_min = obj_diff;
01084 }
01085 }
01086 }
01087 }
01088 else
01089 {
01090 if (!is_upper_bound(j))
01091 {
01092 double grad_diff=Gmaxn-G[j];
01093 if (-G[j] >= Gmaxn2)
01094 Gmaxn2 = -G[j];
01095 if (grad_diff > 0)
01096 {
01097 double obj_diff;
01098 double quad_coef = QD[in]+QD[j]-2*Q_in[j];
01099 if (quad_coef > 0)
01100 obj_diff = -(grad_diff*grad_diff)/quad_coef;
01101 else
01102 obj_diff = -(grad_diff*grad_diff)/TAU;
01103
01104 if (obj_diff <= obj_diff_min)
01105 {
01106 Gmin_idx=j;
01107 obj_diff_min = obj_diff;
01108 }
01109 }
01110 }
01111 }
01112 }
01113
01114 if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
01115 return 1;
01116
01117 if (y[Gmin_idx] == +1)
01118 out_i = Gmaxp_idx;
01119 else
01120 out_i = Gmaxn_idx;
01121 out_j = Gmin_idx;
01122
01123 return 0;
01124 }
01125
01126 bool Solver_NU::be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)
01127 {
01128 if(is_upper_bound(i))
01129 {
01130 if(y[i]==+1)
01131 return(-G[i] > Gmax1);
01132 else
01133 return(-G[i] > Gmax4);
01134 }
01135 else if(is_lower_bound(i))
01136 {
01137 if(y[i]==+1)
01138 return(G[i] > Gmax2);
01139 else
01140 return(G[i] > Gmax3);
01141 }
01142 else
01143 return(false);
01144 }
01145
01146 void Solver_NU::do_shrinking()
01147 {
01148 double Gmax1 = -INF;
01149 double Gmax2 = -INF;
01150 double Gmax3 = -INF;
01151 double Gmax4 = -INF;
01152
01153
01154 int i;
01155 for(i=0;i<active_size;i++)
01156 {
01157 if(!is_upper_bound(i))
01158 {
01159 if(y[i]==+1)
01160 {
01161 if(-G[i] > Gmax1) Gmax1 = -G[i];
01162 }
01163 else if(-G[i] > Gmax4) Gmax4 = -G[i];
01164 }
01165 if(!is_lower_bound(i))
01166 {
01167 if(y[i]==+1)
01168 {
01169 if(G[i] > Gmax2) Gmax2 = G[i];
01170 }
01171 else if(G[i] > Gmax3) Gmax3 = G[i];
01172 }
01173 }
01174
01175 if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10)
01176 {
01177 unshrink = true;
01178 reconstruct_gradient();
01179 active_size = l;
01180 }
01181
01182 for(i=0;i<active_size;i++)
01183 if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
01184 {
01185 active_size--;
01186 while (active_size > i)
01187 {
01188 if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
01189 {
01190 swap_index(i,active_size);
01191 break;
01192 }
01193 active_size--;
01194 }
01195 }
01196 }
01197
01198 double Solver_NU::calculate_rho()
01199 {
01200 int nr_free1 = 0,nr_free2 = 0;
01201 double ub1 = INF, ub2 = INF;
01202 double lb1 = -INF, lb2 = -INF;
01203 double sum_free1 = 0, sum_free2 = 0;
01204
01205 for(int i=0;i<active_size;i++)
01206 {
01207 if(y[i]==+1)
01208 {
01209 if(is_upper_bound(i))
01210 lb1 = max(lb1,G[i]);
01211 else if(is_lower_bound(i))
01212 ub1 = min(ub1,G[i]);
01213 else
01214 {
01215 ++nr_free1;
01216 sum_free1 += G[i];
01217 }
01218 }
01219 else
01220 {
01221 if(is_upper_bound(i))
01222 lb2 = max(lb2,G[i]);
01223 else if(is_lower_bound(i))
01224 ub2 = min(ub2,G[i]);
01225 else
01226 {
01227 ++nr_free2;
01228 sum_free2 += G[i];
01229 }
01230 }
01231 }
01232
01233 double r1,r2;
01234 if(nr_free1 > 0)
01235 r1 = sum_free1/nr_free1;
01236 else
01237 r1 = (ub1+lb1)/2;
01238
01239 if(nr_free2 > 0)
01240 r2 = sum_free2/nr_free2;
01241 else
01242 r2 = (ub2+lb2)/2;
01243
01244 si->r = (r1+r2)/2;
01245 return (r1-r2)/2;
01246 }
01247
01248
01249
01250
01251 class SVC_Q: public Kernel
01252 {
01253 public:
01254 SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
01255 :Kernel(prob.l, prob.x, param)
01256 {
01257 clone(y,y_,prob.l);
01258 cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
01259 QD = new double[prob.l];
01260 for(int i=0;i<prob.l;i++)
01261 QD[i] = (this->*kernel_function)(i,i);
01262 }
01263
01264 Qfloat *get_Q(int i, int len) const
01265 {
01266 Qfloat *data;
01267 int start, j;
01268 if((start = cache->get_data(i,&data,len)) < len)
01269 {
01270 for(j=start;j<len;j++)
01271 data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
01272 }
01273 return data;
01274 }
01275
01276 double *get_QD() const
01277 {
01278 return QD;
01279 }
01280
01281 void swap_index(int i, int j) const
01282 {
01283 cache->swap_index(i,j);
01284 Kernel::swap_index(i,j);
01285 swap(y[i],y[j]);
01286 swap(QD[i],QD[j]);
01287 }
01288
01289 ~SVC_Q()
01290 {
01291 delete[] y;
01292 delete cache;
01293 delete[] QD;
01294 }
01295 private:
01296 schar *y;
01297 Cache *cache;
01298 double *QD;
01299 };
01300
01301 class ONE_CLASS_Q: public Kernel
01302 {
01303 public:
01304 ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
01305 :Kernel(prob.l, prob.x, param)
01306 {
01307 cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
01308 QD = new double[prob.l];
01309 for(int i=0;i<prob.l;i++)
01310 QD[i] = (this->*kernel_function)(i,i);
01311 }
01312
01313 Qfloat *get_Q(int i, int len) const
01314 {
01315 Qfloat *data;
01316 int start, j;
01317 if((start = cache->get_data(i,&data,len)) < len)
01318 {
01319 for(j=start;j<len;j++)
01320 data[j] = (Qfloat)(this->*kernel_function)(i,j);
01321 }
01322 return data;
01323 }
01324
01325 double *get_QD() const
01326 {
01327 return QD;
01328 }
01329
01330 void swap_index(int i, int j) const
01331 {
01332 cache->swap_index(i,j);
01333 Kernel::swap_index(i,j);
01334 swap(QD[i],QD[j]);
01335 }
01336
01337 ~ONE_CLASS_Q()
01338 {
01339 delete cache;
01340 delete[] QD;
01341 }
01342 private:
01343 Cache *cache;
01344 double *QD;
01345 };
01346
01347 class SVR_Q: public Kernel
01348 {
01349 public:
01350 SVR_Q(const svm_problem& prob, const svm_parameter& param)
01351 :Kernel(prob.l, prob.x, param)
01352 {
01353 l = prob.l;
01354 cache = new Cache(l,(long int)(param.cache_size*(1<<20)));
01355 QD = new double[2*l];
01356 sign = new schar[2*l];
01357 index = new int[2*l];
01358 for(int k=0;k<l;k++)
01359 {
01360 sign[k] = 1;
01361 sign[k+l] = -1;
01362 index[k] = k;
01363 index[k+l] = k;
01364 QD[k] = (this->*kernel_function)(k,k);
01365 QD[k+l] = QD[k];
01366 }
01367 buffer[0] = new Qfloat[2*l];
01368 buffer[1] = new Qfloat[2*l];
01369 next_buffer = 0;
01370 }
01371
01372 void swap_index(int i, int j) const
01373 {
01374 swap(sign[i],sign[j]);
01375 swap(index[i],index[j]);
01376 swap(QD[i],QD[j]);
01377 }
01378
01379 Qfloat *get_Q(int i, int len) const
01380 {
01381 Qfloat *data;
01382 int j, real_i = index[i];
01383 if(cache->get_data(real_i,&data,l) < l)
01384 {
01385 for(j=0;j<l;j++)
01386 data[j] = (Qfloat)(this->*kernel_function)(real_i,j);
01387 }
01388
01389
01390 Qfloat *buf = buffer[next_buffer];
01391 next_buffer = 1 - next_buffer;
01392 schar si = sign[i];
01393 for(j=0;j<len;j++)
01394 buf[j] = (Qfloat) si * (Qfloat) sign[j] * data[index[j]];
01395 return buf;
01396 }
01397
01398 double *get_QD() const
01399 {
01400 return QD;
01401 }
01402
01403 ~SVR_Q()
01404 {
01405 delete cache;
01406 delete[] sign;
01407 delete[] index;
01408 delete[] buffer[0];
01409 delete[] buffer[1];
01410 delete[] QD;
01411 }
01412 private:
01413 int l;
01414 Cache *cache;
01415 schar *sign;
01416 int *index;
01417 mutable int next_buffer;
01418 Qfloat *buffer[2];
01419 double *QD;
01420 };
01421
01422
01423
01424
01425 static void solve_c_svc(
01426 const svm_problem *prob, const svm_parameter* param,
01427 double *alpha, Solver::SolutionInfo* si, double Cp, double Cn)
01428 {
01429 int l = prob->l;
01430 double *minus_ones = new double[l];
01431 schar *y = new schar[l];
01432
01433 int i;
01434
01435 for(i=0;i<l;i++)
01436 {
01437 alpha[i] = 0;
01438 minus_ones[i] = -1;
01439 if(prob->y[i] > 0) y[i] = +1; else y[i] = -1;
01440 }
01441
01442 Solver s;
01443 s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
01444 alpha, Cp, Cn, param->eps, si, param->shrinking);
01445
01446 double sum_alpha=0;
01447 for(i=0;i<l;i++)
01448 sum_alpha += alpha[i];
01449
01450 if (Cp==Cn)
01451 info("nu = %f\n", sum_alpha/(Cp*prob->l));
01452
01453 for(i=0;i<l;i++)
01454 alpha[i] *= y[i];
01455
01456 delete[] minus_ones;
01457 delete[] y;
01458 }
01459
01460 static void solve_nu_svc(
01461 const svm_problem *prob, const svm_parameter *param,
01462 double *alpha, Solver::SolutionInfo* si)
01463 {
01464 int i;
01465 int l = prob->l;
01466 double nu = param->nu;
01467
01468 schar *y = new schar[l];
01469
01470 for(i=0;i<l;i++)
01471 if(prob->y[i]>0)
01472 y[i] = +1;
01473 else
01474 y[i] = -1;
01475
01476 double sum_pos = nu*l/2;
01477 double sum_neg = nu*l/2;
01478
01479 for(i=0;i<l;i++)
01480 if(y[i] == +1)
01481 {
01482 alpha[i] = min(1.0,sum_pos);
01483 sum_pos -= alpha[i];
01484 }
01485 else
01486 {
01487 alpha[i] = min(1.0,sum_neg);
01488 sum_neg -= alpha[i];
01489 }
01490
01491 double *zeros = new double[l];
01492
01493 for(i=0;i<l;i++)
01494 zeros[i] = 0;
01495
01496 Solver_NU s;
01497 s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
01498 alpha, 1.0, 1.0, param->eps, si, param->shrinking);
01499 double r = si->r;
01500
01501 info("C = %f\n",1/r);
01502
01503 for(i=0;i<l;i++)
01504 alpha[i] *= y[i]/r;
01505
01506 si->rho /= r;
01507 si->obj /= (r*r);
01508 si->upper_bound_p = 1/r;
01509 si->upper_bound_n = 1/r;
01510
01511 delete[] y;
01512 delete[] zeros;
01513 }
01514
01515 static void solve_one_class(
01516 const svm_problem *prob, const svm_parameter *param,
01517 double *alpha, Solver::SolutionInfo* si)
01518 {
01519 int l = prob->l;
01520 double *zeros = new double[l];
01521 schar *ones = new schar[l];
01522 int i;
01523
01524 int n = (int)(param->nu*prob->l);
01525
01526 for(i=0;i<n;i++)
01527 alpha[i] = 1;
01528 if(n<prob->l)
01529 alpha[n] = param->nu * prob->l - n;
01530 for(i=n+1;i<l;i++)
01531 alpha[i] = 0;
01532
01533 for(i=0;i<l;i++)
01534 {
01535 zeros[i] = 0;
01536 ones[i] = 1;
01537 }
01538
01539 Solver s;
01540 s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
01541 alpha, 1.0, 1.0, param->eps, si, param->shrinking);
01542
01543 delete[] zeros;
01544 delete[] ones;
01545 }
01546
01547 static void solve_epsilon_svr(
01548 const svm_problem *prob, const svm_parameter *param,
01549 double *alpha, Solver::SolutionInfo* si)
01550 {
01551 int l = prob->l;
01552 double *alpha2 = new double[2*l];
01553 double *linear_term = new double[2*l];
01554 schar *y = new schar[2*l];
01555 int i;
01556
01557 for(i=0;i<l;i++)
01558 {
01559 alpha2[i] = 0;
01560 linear_term[i] = param->p - prob->y[i];
01561 y[i] = 1;
01562
01563 alpha2[i+l] = 0;
01564 linear_term[i+l] = param->p + prob->y[i];
01565 y[i+l] = -1;
01566 }
01567
01568 Solver s;
01569 s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
01570 alpha2, param->C, param->C, param->eps, si, param->shrinking);
01571
01572 double sum_alpha = 0;
01573 for(i=0;i<l;i++)
01574 {
01575 alpha[i] = alpha2[i] - alpha2[i+l];
01576 sum_alpha += fabs(alpha[i]);
01577 }
01578 info("nu = %f\n",sum_alpha/(param->C*l));
01579
01580 delete[] alpha2;
01581 delete[] linear_term;
01582 delete[] y;
01583 }
01584
01585 static void solve_nu_svr(
01586 const svm_problem *prob, const svm_parameter *param,
01587 double *alpha, Solver::SolutionInfo* si)
01588 {
01589 int l = prob->l;
01590 double C = param->C;
01591 double *alpha2 = new double[2*l];
01592 double *linear_term = new double[2*l];
01593 schar *y = new schar[2*l];
01594 int i;
01595
01596 double sum = C * param->nu * l / 2;
01597 for(i=0;i<l;i++)
01598 {
01599 alpha2[i] = alpha2[i+l] = min(sum,C);
01600 sum -= alpha2[i];
01601
01602 linear_term[i] = - prob->y[i];
01603 y[i] = 1;
01604
01605 linear_term[i+l] = prob->y[i];
01606 y[i+l] = -1;
01607 }
01608
01609 Solver_NU s;
01610 s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
01611 alpha2, C, C, param->eps, si, param->shrinking);
01612
01613 info("epsilon = %f\n",-si->r);
01614
01615 for(i=0;i<l;i++)
01616 alpha[i] = alpha2[i] - alpha2[i+l];
01617
01618 delete[] alpha2;
01619 delete[] linear_term;
01620 delete[] y;
01621 }
01622
01623
01624
01625
01626 struct decision_function
01627 {
01628 double *alpha;
01629 double rho;
01630 };
01631
01632 static decision_function svm_train_one(
01633 const svm_problem *prob, const svm_parameter *param,
01634 double Cp, double Cn)
01635 {
01636 double *alpha = Malloc(double,prob->l);
01637 Solver::SolutionInfo si;
01638 switch(param->svm_type)
01639 {
01640 case C_SVC:
01641 solve_c_svc(prob,param,alpha,&si,Cp,Cn);
01642 break;
01643 case NU_SVC:
01644 solve_nu_svc(prob,param,alpha,&si);
01645 break;
01646 case ONE_CLASS:
01647 solve_one_class(prob,param,alpha,&si);
01648 break;
01649 case EPSILON_SVR:
01650 solve_epsilon_svr(prob,param,alpha,&si);
01651 break;
01652 case NU_SVR:
01653 solve_nu_svr(prob,param,alpha,&si);
01654 break;
01655 }
01656
01657 info("obj = %f, rho = %f\n",si.obj,si.rho);
01658
01659
01660
01661 int nSV = 0;
01662 int nBSV = 0;
01663 for(int i=0;i<prob->l;i++)
01664 {
01665 if(fabs(alpha[i]) > 0)
01666 {
01667 ++nSV;
01668 if(prob->y[i] > 0)
01669 {
01670 if(fabs(alpha[i]) >= si.upper_bound_p)
01671 ++nBSV;
01672 }
01673 else
01674 {
01675 if(fabs(alpha[i]) >= si.upper_bound_n)
01676 ++nBSV;
01677 }
01678 }
01679 }
01680
01681 info("nSV = %d, nBSV = %d\n",nSV,nBSV);
01682
01683 decision_function f;
01684 f.alpha = alpha;
01685 f.rho = si.rho;
01686 return f;
01687 }
01688
01689
01690 static void sigmoid_train(
01691 int l, const double *dec_values, const double *labels,
01692 double& A, double& B)
01693 {
01694 double prior1=0, prior0 = 0;
01695 int i;
01696
01697 for (i=0;i<l;i++)
01698 if (labels[i] > 0) prior1+=1;
01699 else prior0+=1;
01700
01701 int max_iter=100;
01702 double min_step=1e-10;
01703 double sigma=1e-12;
01704 double eps=1e-5;
01705 double hiTarget=(prior1+1.0)/(prior1+2.0);
01706 double loTarget=1/(prior0+2.0);
01707 double *t=Malloc(double,l);
01708 double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
01709 double newA,newB,newf,d1,d2;
01710 int iter;
01711
01712
01713 A=0.0; B=log((prior0+1.0)/(prior1+1.0));
01714 double fval = 0.0;
01715
01716 for (i=0;i<l;i++)
01717 {
01718 if (labels[i]>0) t[i]=hiTarget;
01719 else t[i]=loTarget;
01720 fApB = dec_values[i]*A+B;
01721 if (fApB>=0)
01722 fval += t[i]*fApB + log(1+exp(-fApB));
01723 else
01724 fval += (t[i] - 1)*fApB +log(1+exp(fApB));
01725 }
01726 for (iter=0;iter<max_iter;iter++)
01727 {
01728
01729 h11=sigma;
01730 h22=sigma;
01731 h21=0.0;g1=0.0;g2=0.0;
01732 for (i=0;i<l;i++)
01733 {
01734 fApB = dec_values[i]*A+B;
01735 if (fApB >= 0)
01736 {
01737 p=exp(-fApB)/(1.0+exp(-fApB));
01738 q=1.0/(1.0+exp(-fApB));
01739 }
01740 else
01741 {
01742 p=1.0/(1.0+exp(fApB));
01743 q=exp(fApB)/(1.0+exp(fApB));
01744 }
01745 d2=p*q;
01746 h11+=dec_values[i]*dec_values[i]*d2;
01747 h22+=d2;
01748 h21+=dec_values[i]*d2;
01749 d1=t[i]-p;
01750 g1+=dec_values[i]*d1;
01751 g2+=d1;
01752 }
01753
01754
01755 if (fabs(g1)<eps && fabs(g2)<eps)
01756 break;
01757
01758
01759 det=h11*h22-h21*h21;
01760 dA=-(h22*g1 - h21 * g2) / det;
01761 dB=-(-h21*g1+ h11 * g2) / det;
01762 gd=g1*dA+g2*dB;
01763
01764
01765 stepsize = 1;
01766 while (stepsize >= min_step)
01767 {
01768 newA = A + stepsize * dA;
01769 newB = B + stepsize * dB;
01770
01771
01772 newf = 0.0;
01773 for (i=0;i<l;i++)
01774 {
01775 fApB = dec_values[i]*newA+newB;
01776 if (fApB >= 0)
01777 newf += t[i]*fApB + log(1+exp(-fApB));
01778 else
01779 newf += (t[i] - 1)*fApB +log(1+exp(fApB));
01780 }
01781
01782 if (newf<fval+0.0001*stepsize*gd)
01783 {
01784 A=newA;B=newB;fval=newf;
01785 break;
01786 }
01787 else
01788 stepsize = stepsize / 2.0;
01789 }
01790
01791 if (stepsize < min_step)
01792 {
01793 info("Line search fails in two-class probability estimates\n");
01794 break;
01795 }
01796 }
01797
01798 if (iter>=max_iter)
01799 info("Reaching maximal iterations in two-class probability estimates\n");
01800 free(t);
01801 }
01802
01803 static double sigmoid_predict(double decision_value, double A, double B)
01804 {
01805 double fApB = decision_value*A+B;
01806 if (fApB >= 0)
01807 return exp(-fApB)/(1.0+exp(-fApB));
01808 else
01809 return 1.0/(1+exp(fApB)) ;
01810 }
01811
01812
01813 static void multiclass_probability(int k, double **r, double *p)
01814 {
01815 int t,j;
01816 int iter = 0, max_iter=max(100,k);
01817 double **Q=Malloc(double *,k);
01818 double *Qp=Malloc(double,k);
01819 double pQp, eps=0.005/k;
01820
01821 for (t=0;t<k;t++)
01822 {
01823 p[t]=1.0/k;
01824 Q[t]=Malloc(double,k);
01825 Q[t][t]=0;
01826 for (j=0;j<t;j++)
01827 {
01828 Q[t][t]+=r[j][t]*r[j][t];
01829 Q[t][j]=Q[j][t];
01830 }
01831 for (j=t+1;j<k;j++)
01832 {
01833 Q[t][t]+=r[j][t]*r[j][t];
01834 Q[t][j]=-r[j][t]*r[t][j];
01835 }
01836 }
01837 for (iter=0;iter<max_iter;iter++)
01838 {
01839
01840 pQp=0;
01841 for (t=0;t<k;t++)
01842 {
01843 Qp[t]=0;
01844 for (j=0;j<k;j++)
01845 Qp[t]+=Q[t][j]*p[j];
01846 pQp+=p[t]*Qp[t];
01847 }
01848 double max_error=0;
01849 for (t=0;t<k;t++)
01850 {
01851 double error=fabs(Qp[t]-pQp);
01852 if (error>max_error)
01853 max_error=error;
01854 }
01855 if (max_error<eps) break;
01856
01857 for (t=0;t<k;t++)
01858 {
01859 double diff=(-Qp[t]+pQp)/Q[t][t];
01860 p[t]+=diff;
01861 pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
01862 for (j=0;j<k;j++)
01863 {
01864 Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
01865 p[j]/=(1+diff);
01866 }
01867 }
01868 }
01869 if (iter>=max_iter)
01870 info("Exceeds max_iter in multiclass_prob\n");
01871 for(t=0;t<k;t++) free(Q[t]);
01872 free(Q);
01873 free(Qp);
01874 }
01875
01876
01877 static void svm_binary_svc_probability(
01878 const svm_problem *prob, const svm_parameter *param,
01879 double Cp, double Cn, double& probA, double& probB)
01880 {
01881 int i;
01882 int nr_fold = 5;
01883 int *perm = Malloc(int,prob->l);
01884 double *dec_values = Malloc(double,prob->l);
01885
01886
01887 for(i=0;i<prob->l;i++) perm[i]=i;
01888 for(i=0;i<prob->l;i++)
01889 {
01890 int j = i+rand()%(prob->l-i);
01891 swap(perm[i],perm[j]);
01892 }
01893 for(i=0;i<nr_fold;i++)
01894 {
01895 int begin = i*prob->l/nr_fold;
01896 int end = (i+1)*prob->l/nr_fold;
01897 int j,k;
01898 struct svm_problem subprob;
01899
01900 subprob.l = prob->l-(end-begin);
01901 subprob.x = Malloc(struct svm_node*,subprob.l);
01902 subprob.y = Malloc(double,subprob.l);
01903
01904 k=0;
01905 for(j=0;j<begin;j++)
01906 {
01907 subprob.x[k] = prob->x[perm[j]];
01908 subprob.y[k] = prob->y[perm[j]];
01909 ++k;
01910 }
01911 for(j=end;j<prob->l;j++)
01912 {
01913 subprob.x[k] = prob->x[perm[j]];
01914 subprob.y[k] = prob->y[perm[j]];
01915 ++k;
01916 }
01917 int p_count=0,n_count=0;
01918 for(j=0;j<k;j++)
01919 if(subprob.y[j]>0)
01920 p_count++;
01921 else
01922 n_count++;
01923
01924 if(p_count==0 && n_count==0)
01925 for(j=begin;j<end;j++)
01926 dec_values[perm[j]] = 0;
01927 else if(p_count > 0 && n_count == 0)
01928 for(j=begin;j<end;j++)
01929 dec_values[perm[j]] = 1;
01930 else if(p_count == 0 && n_count > 0)
01931 for(j=begin;j<end;j++)
01932 dec_values[perm[j]] = -1;
01933 else
01934 {
01935 svm_parameter subparam = *param;
01936 subparam.probability=0;
01937 subparam.C=1.0;
01938 subparam.nr_weight=2;
01939 subparam.weight_label = Malloc(int,2);
01940 subparam.weight = Malloc(double,2);
01941 subparam.weight_label[0]=+1;
01942 subparam.weight_label[1]=-1;
01943 subparam.weight[0]=Cp;
01944 subparam.weight[1]=Cn;
01945 struct svm_model *submodel = svm_train(&subprob,&subparam);
01946 for(j=begin;j<end;j++)
01947 {
01948 svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]]));
01949
01950 dec_values[perm[j]] *= submodel->label[0];
01951 }
01952 svm_free_and_destroy_model(&submodel);
01953 svm_destroy_param(&subparam);
01954 }
01955 free(subprob.x);
01956 free(subprob.y);
01957 }
01958 sigmoid_train(prob->l,dec_values,prob->y,probA,probB);
01959 free(dec_values);
01960 free(perm);
01961 }
01962
01963
01964 static double svm_svr_probability(
01965 const svm_problem *prob, const svm_parameter *param)
01966 {
01967 int i;
01968 int nr_fold = 5;
01969 double *ymv = Malloc(double,prob->l);
01970 double mae = 0;
01971
01972 svm_parameter newparam = *param;
01973 newparam.probability = 0;
01974 svm_cross_validation(prob,&newparam,nr_fold,ymv);
01975 for(i=0;i<prob->l;i++)
01976 {
01977 ymv[i]=prob->y[i]-ymv[i];
01978 mae += fabs(ymv[i]);
01979 }
01980 mae /= prob->l;
01981 double std=sqrt(2*mae*mae);
01982 int count=0;
01983 mae=0;
01984 for(i=0;i<prob->l;i++)
01985 if (fabs(ymv[i]) > 5*std)
01986 count=count+1;
01987 else
01988 mae+=fabs(ymv[i]);
01989 mae /= (prob->l-count);
01990 info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
01991 free(ymv);
01992 return mae;
01993 }
01994
01995
01996
01997
01998 static void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)
01999 {
02000 int l = prob->l;
02001 int max_nr_class = 16;
02002 int nr_class = 0;
02003 int *label = Malloc(int,max_nr_class);
02004 int *count = Malloc(int,max_nr_class);
02005 int *data_label = Malloc(int,l);
02006 int i;
02007
02008 for(i=0;i<l;i++)
02009 {
02010 int this_label = (int)prob->y[i];
02011 int j;
02012 for(j=0;j<nr_class;j++)
02013 {
02014 if(this_label == label[j])
02015 {
02016 ++count[j];
02017 break;
02018 }
02019 }
02020 data_label[i] = j;
02021 if(j == nr_class)
02022 {
02023 if(nr_class == max_nr_class)
02024 {
02025 max_nr_class *= 2;
02026 label = (int *)realloc(label,max_nr_class*sizeof(int));
02027 count = (int *)realloc(count,max_nr_class*sizeof(int));
02028 }
02029 label[nr_class] = this_label;
02030 count[nr_class] = 1;
02031 ++nr_class;
02032 }
02033 }
02034
02035 int *start = Malloc(int,nr_class);
02036 start[0] = 0;
02037 for(i=1;i<nr_class;i++)
02038 start[i] = start[i-1]+count[i-1];
02039 for(i=0;i<l;i++)
02040 {
02041 perm[start[data_label[i]]] = i;
02042 ++start[data_label[i]];
02043 }
02044 start[0] = 0;
02045 for(i=1;i<nr_class;i++)
02046 start[i] = start[i-1]+count[i-1];
02047
02048 *nr_class_ret = nr_class;
02049 *label_ret = label;
02050 *start_ret = start;
02051 *count_ret = count;
02052 free(data_label);
02053 }
02054
02055
02056
02057
02058 svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
02059 {
02060 svm_model *model = Malloc(svm_model,1);
02061 model->param = *param;
02062 model->free_sv = 0;
02063
02064 if(param->svm_type == ONE_CLASS ||
02065 param->svm_type == EPSILON_SVR ||
02066 param->svm_type == NU_SVR)
02067 {
02068
02069 model->nr_class = 2;
02070 model->label = NULL;
02071 model->nSV = NULL;
02072 model->probA = NULL; model->probB = NULL;
02073 model->sv_coef = Malloc(double *,1);
02074
02075 if(param->probability &&
02076 (param->svm_type == EPSILON_SVR ||
02077 param->svm_type == NU_SVR))
02078 {
02079 model->probA = Malloc(double,1);
02080 model->probA[0] = svm_svr_probability(prob,param);
02081 }
02082
02083 decision_function f = svm_train_one(prob,param,0,0);
02084 model->rho = Malloc(double,1);
02085 model->rho[0] = f.rho;
02086
02087 int nSV = 0;
02088 int i;
02089 for(i=0;i<prob->l;i++)
02090 if(fabs(f.alpha[i]) > 0) ++nSV;
02091 model->l = nSV;
02092 model->SV = Malloc(svm_node *,nSV);
02093 model->sv_coef[0] = Malloc(double,nSV);
02094 int j = 0;
02095 for(i=0;i<prob->l;i++)
02096 if(fabs(f.alpha[i]) > 0)
02097 {
02098 model->SV[j] = prob->x[i];
02099 model->sv_coef[0][j] = f.alpha[i];
02100 ++j;
02101 }
02102
02103 free(f.alpha);
02104 }
02105 else
02106 {
02107
02108 int l = prob->l;
02109 int nr_class;
02110 int *label = NULL;
02111 int *start = NULL;
02112 int *count = NULL;
02113 int *perm = Malloc(int,l);
02114
02115
02116 svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
02117 svm_node **x = Malloc(svm_node *,l);
02118 int i;
02119 for(i=0;i<l;i++)
02120 x[i] = prob->x[perm[i]];
02121
02122
02123
02124 double *weighted_C = Malloc(double, nr_class);
02125 for(i=0;i<nr_class;i++)
02126 weighted_C[i] = param->C;
02127 for(i=0;i<param->nr_weight;i++)
02128 {
02129 int j;
02130 for(j=0;j<nr_class;j++)
02131 if(param->weight_label[i] == label[j])
02132 break;
02133 if(j == nr_class)
02134 fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]);
02135 else
02136 weighted_C[j] *= param->weight[i];
02137 }
02138
02139
02140
02141 bool *nonzero = Malloc(bool,l);
02142 for(i=0;i<l;i++)
02143 nonzero[i] = false;
02144 decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
02145
02146 double *probA=NULL,*probB=NULL;
02147 if (param->probability)
02148 {
02149 probA=Malloc(double,nr_class*(nr_class-1)/2);
02150 probB=Malloc(double,nr_class*(nr_class-1)/2);
02151 }
02152
02153 int p = 0;
02154 for(i=0;i<nr_class;i++)
02155 for(int j=i+1;j<nr_class;j++)
02156 {
02157 svm_problem sub_prob;
02158 int si = start[i], sj = start[j];
02159 int ci = count[i], cj = count[j];
02160 sub_prob.l = ci+cj;
02161 sub_prob.x = Malloc(svm_node *,sub_prob.l);
02162 sub_prob.y = Malloc(double,sub_prob.l);
02163 int k;
02164 for(k=0;k<ci;k++)
02165 {
02166 sub_prob.x[k] = x[si+k];
02167 sub_prob.y[k] = +1;
02168 }
02169 for(k=0;k<cj;k++)
02170 {
02171 sub_prob.x[ci+k] = x[sj+k];
02172 sub_prob.y[ci+k] = -1;
02173 }
02174
02175 if(param->probability)
02176 svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);
02177
02178 f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
02179 for(k=0;k<ci;k++)
02180 if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
02181 nonzero[si+k] = true;
02182 for(k=0;k<cj;k++)
02183 if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
02184 nonzero[sj+k] = true;
02185 free(sub_prob.x);
02186 free(sub_prob.y);
02187 ++p;
02188 }
02189
02190
02191
02192 model->nr_class = nr_class;
02193
02194 model->label = Malloc(int,nr_class);
02195 for(i=0;i<nr_class;i++)
02196 model->label[i] = label[i];
02197
02198 model->rho = Malloc(double,nr_class*(nr_class-1)/2);
02199 for(i=0;i<nr_class*(nr_class-1)/2;i++)
02200 model->rho[i] = f[i].rho;
02201
02202 if(param->probability)
02203 {
02204 model->probA = Malloc(double,nr_class*(nr_class-1)/2);
02205 model->probB = Malloc(double,nr_class*(nr_class-1)/2);
02206 for(i=0;i<nr_class*(nr_class-1)/2;i++)
02207 {
02208 model->probA[i] = probA[i];
02209 model->probB[i] = probB[i];
02210 }
02211 }
02212 else
02213 {
02214 model->probA=NULL;
02215 model->probB=NULL;
02216 }
02217
02218 int total_sv = 0;
02219 int *nz_count = Malloc(int,nr_class);
02220 model->nSV = Malloc(int,nr_class);
02221 for(i=0;i<nr_class;i++)
02222 {
02223 int nSV = 0;
02224 for(int j=0;j<count[i];j++)
02225 if(nonzero[start[i]+j])
02226 {
02227 ++nSV;
02228 ++total_sv;
02229 }
02230 model->nSV[i] = nSV;
02231 nz_count[i] = nSV;
02232 }
02233
02234 info("Total nSV = %d\n",total_sv);
02235
02236 model->l = total_sv;
02237 model->SV = Malloc(svm_node *,total_sv);
02238 p = 0;
02239 for(i=0;i<l;i++)
02240 if(nonzero[i])
02241 {
02242 model->SV[p++] = x[i];
02243
02244 }
02245
02246
02247 int *nz_start = Malloc(int,nr_class);
02248 nz_start[0] = 0;
02249 for(i=1;i<nr_class;i++)
02250 nz_start[i] = nz_start[i-1]+nz_count[i-1];
02251
02252 model->sv_coef = Malloc(double *,nr_class-1);
02253 for(i=0;i<nr_class-1;i++)
02254 model->sv_coef[i] = Malloc(double,total_sv);
02255
02256 p = 0;
02257 for(i=0;i<nr_class;i++)
02258 for(int j=i+1;j<nr_class;j++)
02259 {
02260
02261
02262
02263
02264 int si = start[i];
02265 int sj = start[j];
02266 int ci = count[i];
02267 int cj = count[j];
02268
02269 int q = nz_start[i];
02270 int k;
02271 for(k=0;k<ci;k++)
02272 if(nonzero[si+k])
02273 model->sv_coef[j-1][q++] = f[p].alpha[k];
02274 q = nz_start[j];
02275 for(k=0;k<cj;k++)
02276 if(nonzero[sj+k])
02277 model->sv_coef[i][q++] = f[p].alpha[ci+k];
02278 ++p;
02279 }
02280
02281 model->SVidx = (unsigned char*) nonzero;
02282 free(label);
02283 free(probA);
02284 free(probB);
02285 free(count);
02286 free(perm);
02287 free(start);
02288 free(x);
02289 free(weighted_C);
02290
02291 for(i=0;i<nr_class*(nr_class-1)/2;i++)
02292 free(f[i].alpha);
02293 free(f);
02294 free(nz_count);
02295 free(nz_start);
02296 }
02297 return model;
02298 }
02299
02300
02301 void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target)
02302 {
02303 int i;
02304 int *fold_start = Malloc(int,nr_fold+1);
02305 int l = prob->l;
02306 int *perm = Malloc(int,l);
02307 int nr_class;
02308
02309
02310
02311 if((param->svm_type == C_SVC ||
02312 param->svm_type == NU_SVC) && nr_fold < l)
02313 {
02314 int *start = NULL;
02315 int *label = NULL;
02316 int *count = NULL;
02317 svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
02318
02319
02320 int *fold_count = Malloc(int,nr_fold);
02321 int c;
02322 int *index = Malloc(int,l);
02323 for(i=0;i<l;i++)
02324 index[i]=perm[i];
02325 for (c=0; c<nr_class; c++)
02326 for(i=0;i<count[c];i++)
02327 {
02328 int j = i+rand()%(count[c]-i);
02329 swap(index[start[c]+j],index[start[c]+i]);
02330 }
02331 for(i=0;i<nr_fold;i++)
02332 {
02333 fold_count[i] = 0;
02334 for (c=0; c<nr_class;c++)
02335 fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold;
02336 }
02337 fold_start[0]=0;
02338 for (i=1;i<=nr_fold;i++)
02339 fold_start[i] = fold_start[i-1]+fold_count[i-1];
02340 for (c=0; c<nr_class;c++)
02341 for(i=0;i<nr_fold;i++)
02342 {
02343 int begin = start[c]+i*count[c]/nr_fold;
02344 int end = start[c]+(i+1)*count[c]/nr_fold;
02345 for(int j=begin;j<end;j++)
02346 {
02347 perm[fold_start[i]] = index[j];
02348 fold_start[i]++;
02349 }
02350 }
02351 fold_start[0]=0;
02352 for (i=1;i<=nr_fold;i++)
02353 fold_start[i] = fold_start[i-1]+fold_count[i-1];
02354 free(start);
02355 free(label);
02356 free(count);
02357 free(index);
02358 free(fold_count);
02359 }
02360 else
02361 {
02362 for(i=0;i<l;i++) perm[i]=i;
02363 for(i=0;i<l;i++)
02364 {
02365 int j = i+rand()%(l-i);
02366 swap(perm[i],perm[j]);
02367 }
02368 for(i=0;i<=nr_fold;i++)
02369 fold_start[i]=i*l/nr_fold;
02370 }
02371
02372 for(i=0;i<nr_fold;i++)
02373 {
02374 int begin = fold_start[i];
02375 int end = fold_start[i+1];
02376 int j,k;
02377 struct svm_problem subprob;
02378
02379 subprob.l = l-(end-begin);
02380 subprob.x = Malloc(struct svm_node*,subprob.l);
02381 subprob.y = Malloc(double,subprob.l);
02382
02383 k=0;
02384 for(j=0;j<begin;j++)
02385 {
02386 subprob.x[k] = prob->x[perm[j]];
02387 subprob.y[k] = prob->y[perm[j]];
02388 ++k;
02389 }
02390 for(j=end;j<l;j++)
02391 {
02392 subprob.x[k] = prob->x[perm[j]];
02393 subprob.y[k] = prob->y[perm[j]];
02394 ++k;
02395 }
02396 struct svm_model *submodel = svm_train(&subprob,param);
02397 if(param->probability &&
02398 (param->svm_type == C_SVC || param->svm_type == NU_SVC))
02399 {
02400 double *prob_estimates=Malloc(double,svm_get_nr_class(submodel));
02401 for(j=begin;j<end;j++)
02402 target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
02403 free(prob_estimates);
02404 }
02405 else
02406 for(j=begin;j<end;j++)
02407 target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
02408 svm_free_and_destroy_model(&submodel);
02409 free(subprob.x);
02410 free(subprob.y);
02411 }
02412 free(fold_start);
02413 free(perm);
02414 }
02415
02416
02417 int svm_get_svm_type(const svm_model *model)
02418 {
02419 return model->param.svm_type;
02420 }
02421
02422 int svm_get_nr_class(const svm_model *model)
02423 {
02424 return model->nr_class;
02425 }
02426
02427 void svm_get_labels(const svm_model *model, int* label)
02428 {
02429 if (model->label != NULL)
02430 for(int i=0;i<model->nr_class;i++)
02431 label[i] = model->label[i];
02432 }
02433
02434 void svm_get_support_vectors(const svm_model *model, int nx, int* sv)
02435 {
02436
02437 if (model->label != NULL)
02438 {
02439 for(int i=0;i<nx;i++)
02440 {
02441
02442 sv[i] = model->SVidx[i];
02443 }
02444 }
02445 }
02446
02447 double svm_get_svr_probability(const svm_model *model)
02448 {
02449 if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
02450 model->probA!=NULL)
02451 return model->probA[0];
02452 else
02453 {
02454 fprintf(stderr,"Model doesn't contain information for SVR probability inference\n");
02455 return 0;
02456 }
02457 }
02458
02459 double svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
02460 {
02461 if(model->param.svm_type == ONE_CLASS ||
02462 model->param.svm_type == EPSILON_SVR ||
02463 model->param.svm_type == NU_SVR)
02464 {
02465 double *sv_coef = model->sv_coef[0];
02466 double sum = 0;
02467 for(int i=0;i<model->l;i++)
02468 sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
02469 sum -= model->rho[0];
02470 *dec_values = sum;
02471
02472 if(model->param.svm_type == ONE_CLASS)
02473 return (sum>0)?1:-1;
02474 else
02475 return sum;
02476 }
02477 else
02478 {
02479 int i;
02480 int nr_class = model->nr_class;
02481 int l = model->l;
02482
02483 double *kvalue = Malloc(double,l);
02484 for(i=0;i<l;i++)
02485 kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);
02486
02487 int *start = Malloc(int,nr_class);
02488 start[0] = 0;
02489 for(i=1;i<nr_class;i++)
02490 start[i] = start[i-1]+model->nSV[i-1];
02491
02492 int *vote = Malloc(int,nr_class);
02493 for(i=0;i<nr_class;i++)
02494 vote[i] = 0;
02495
02496 int p=0;
02497
02498 for(i=0;i<nr_class;i++)
02499 for(int j=i+1;j<nr_class;j++)
02500 {
02501 double sum = 0;
02502 int si = start[i];
02503 int sj = start[j];
02504 int ci = model->nSV[i];
02505 int cj = model->nSV[j];
02506
02507 int k;
02508 double *coef1 = model->sv_coef[j-1];
02509 double *coef2 = model->sv_coef[i];
02510 for(k=0;k<ci;k++)
02511 sum += coef1[si+k] * kvalue[si+k];
02512 for(k=0;k<cj;k++)
02513 sum += coef2[sj+k] * kvalue[sj+k];
02514 sum -= model->rho[p];
02515 dec_values[p] = sum;
02516
02517
02518 if(dec_values[p] > 0)
02519 ++vote[i];
02520 else
02521 ++vote[j];
02522 p++;
02523 }
02524
02525 int vote_max_idx = 0;
02526 for(i=1;i<nr_class;i++)
02527 if(vote[i] > vote[vote_max_idx])
02528 vote_max_idx = i;
02529
02530 free(kvalue);
02531 free(start);
02532 free(vote);
02533 return model->label[vote_max_idx];
02534 }
02535 }
02536
02537 double svm_predict(const svm_model *model, const svm_node *x)
02538 {
02539 int nr_class = model->nr_class;
02540 double *dec_values;
02541 if(model->param.svm_type == ONE_CLASS ||
02542 model->param.svm_type == EPSILON_SVR ||
02543 model->param.svm_type == NU_SVR)
02544 dec_values = Malloc(double, 1);
02545 else
02546 dec_values = Malloc(double, nr_class*(nr_class-1)/2);
02547 double pred_result = svm_predict_values(model, x, dec_values);
02548 free(dec_values);
02549 return pred_result;
02550 }
02551
02552 double svm_predict_probability(
02553 const svm_model *model, const svm_node *x, double *prob_estimates)
02554 {
02555 if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
02556 model->probA!=NULL && model->probB!=NULL)
02557 {
02558 int i;
02559 int nr_class = model->nr_class;
02560 double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
02561 svm_predict_values(model, x, dec_values);
02562
02563 double min_prob=1e-7;
02564 double **pairwise_prob=Malloc(double *,nr_class);
02565 for(i=0;i<nr_class;i++)
02566 pairwise_prob[i]=Malloc(double,nr_class);
02567 int k=0;
02568 for(i=0;i<nr_class;i++)
02569 for(int j=i+1;j<nr_class;j++)
02570 {
02571 pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
02572 pairwise_prob[j][i]=1-pairwise_prob[i][j];
02573 k++;
02574 }
02575 multiclass_probability(nr_class,pairwise_prob,prob_estimates);
02576
02577 int prob_max_idx = 0;
02578 for(i=1;i<nr_class;i++)
02579 if(prob_estimates[i] > prob_estimates[prob_max_idx])
02580 prob_max_idx = i;
02581 for(i=0;i<nr_class;i++)
02582 free(pairwise_prob[i]);
02583 free(dec_values);
02584 free(pairwise_prob);
02585 return model->label[prob_max_idx];
02586 }
02587 else
02588 return svm_predict(model, x);
02589 }
02590
02591 static const char *svm_type_table[] =
02592 {
02593 "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
02594 };
02595
02596 static const char *kernel_type_table[]=
02597 {
02598 "linear","polynomial","rbf","sigmoid","precomputed",NULL
02599 };
02600
02601 int svm_save_model(const char *model_file_name, const svm_model *model)
02602 {
02603 FILE *fp = fopen(model_file_name,"w");
02604 if(fp==NULL) return -1;
02605
02606 const svm_parameter& param = model->param;
02607
02608 fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
02609 fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);
02610
02611 if(param.kernel_type == POLY)
02612 fprintf(fp,"degree %d\n", param.degree);
02613
02614 if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
02615 fprintf(fp,"gamma %g\n", param.gamma);
02616
02617 if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
02618 fprintf(fp,"coef0 %g\n", param.coef0);
02619
02620 int nr_class = model->nr_class;
02621 int l = model->l;
02622 fprintf(fp, "nr_class %d\n", nr_class);
02623 fprintf(fp, "total_sv %d\n",l);
02624
02625 {
02626 fprintf(fp, "rho");
02627 for(int i=0;i<nr_class*(nr_class-1)/2;i++)
02628 fprintf(fp," %g",model->rho[i]);
02629 fprintf(fp, "\n");
02630 }
02631
02632 if(model->label)
02633 {
02634 fprintf(fp, "label");
02635 for(int i=0;i<nr_class;i++)
02636 fprintf(fp," %d",model->label[i]);
02637 fprintf(fp, "\n");
02638 }
02639
02640 if(model->probA)
02641 {
02642 fprintf(fp, "probA");
02643 for(int i=0;i<nr_class*(nr_class-1)/2;i++)
02644 fprintf(fp," %g",model->probA[i]);
02645 fprintf(fp, "\n");
02646 }
02647 if(model->probB)
02648 {
02649 fprintf(fp, "probB");
02650 for(int i=0;i<nr_class*(nr_class-1)/2;i++)
02651 fprintf(fp," %g",model->probB[i]);
02652 fprintf(fp, "\n");
02653 }
02654
02655 if(model->nSV)
02656 {
02657 fprintf(fp, "nr_sv");
02658 for(int i=0;i<nr_class;i++)
02659 fprintf(fp," %d",model->nSV[i]);
02660 fprintf(fp, "\n");
02661 }
02662
02663 fprintf(fp, "SV\n");
02664 const double * const *sv_coef = model->sv_coef;
02665 const svm_node * const *SV = model->SV;
02666
02667 for(int i=0;i<l;i++)
02668 {
02669 for(int j=0;j<nr_class-1;j++)
02670 fprintf(fp, "%.16g ",sv_coef[j][i]);
02671
02672 const svm_node *p = SV[i];
02673
02674 if(param.kernel_type == PRECOMPUTED)
02675 fprintf(fp,"0:%d ",(int)(p->value));
02676 else
02677 while(p->index != -1)
02678 {
02679 fprintf(fp,"%d:%.8g ",p->index,p->value);
02680 p++;
02681 }
02682 fprintf(fp, "\n");
02683 }
02684 if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
02685 else return 0;
02686 }
02687
02688 static char *line = NULL;
02689 static int max_line_len;
02690
02691 static char* readline(FILE *input)
02692 {
02693 int len;
02694
02695 if(fgets(line,max_line_len,input) == NULL)
02696 return NULL;
02697
02698 while(strrchr(line,'\n') == NULL)
02699 {
02700 max_line_len *= 2;
02701 line = (char *) realloc(line,max_line_len);
02702 len = (int) strlen(line);
02703 if(fgets(line+len,max_line_len-len,input) == NULL)
02704 break;
02705 }
02706 return line;
02707 }
02708
02709 svm_model *svm_load_model(const char *model_file_name)
02710 {
02711 FILE *fp = fopen(model_file_name,"rb");
02712 if(fp==NULL) return NULL;
02713
02714
02715
02716 svm_model *model = Malloc(svm_model,1);
02717 svm_parameter& param = model->param;
02718 model->rho = NULL;
02719 model->probA = NULL;
02720 model->probB = NULL;
02721 model->label = NULL;
02722 model->nSV = NULL;
02723
02724 char cmd[81];
02725 while(1)
02726 {
02727 fscanf(fp,"%80s",cmd);
02728
02729 if(strcmp(cmd,"svm_type")==0)
02730 {
02731 fscanf(fp,"%80s",cmd);
02732 int i;
02733 for(i=0;svm_type_table[i];i++)
02734 {
02735 if(strcmp(svm_type_table[i],cmd)==0)
02736 {
02737 param.svm_type=i;
02738 break;
02739 }
02740 }
02741 if(svm_type_table[i] == NULL)
02742 {
02743 fprintf(stderr,"unknown svm type.\n");
02744 free(model->rho);
02745 free(model->label);
02746 free(model->nSV);
02747 free(model);
02748 return NULL;
02749 }
02750 }
02751 else if(strcmp(cmd,"kernel_type")==0)
02752 {
02753 fscanf(fp,"%80s",cmd);
02754 int i;
02755 for(i=0;kernel_type_table[i];i++)
02756 {
02757 if(strcmp(kernel_type_table[i],cmd)==0)
02758 {
02759 param.kernel_type=i;
02760 break;
02761 }
02762 }
02763 if(kernel_type_table[i] == NULL)
02764 {
02765 fprintf(stderr,"unknown kernel function.\n");
02766 free(model->rho);
02767 free(model->label);
02768 free(model->nSV);
02769 free(model);
02770 return NULL;
02771 }
02772 }
02773 else if(strcmp(cmd,"degree")==0)
02774 fscanf(fp,"%d",¶m.degree);
02775 else if(strcmp(cmd,"gamma")==0)
02776 fscanf(fp,"%lf",¶m.gamma);
02777 else if(strcmp(cmd,"coef0")==0)
02778 fscanf(fp,"%lf",¶m.coef0);
02779 else if(strcmp(cmd,"nr_class")==0)
02780 fscanf(fp,"%d",&model->nr_class);
02781 else if(strcmp(cmd,"total_sv")==0)
02782 fscanf(fp,"%d",&model->l);
02783 else if(strcmp(cmd,"rho")==0)
02784 {
02785 int n = model->nr_class * (model->nr_class-1)/2;
02786 model->rho = Malloc(double,n);
02787 for(int i=0;i<n;i++)
02788 fscanf(fp,"%lf",&model->rho[i]);
02789 }
02790 else if(strcmp(cmd,"label")==0)
02791 {
02792 int n = model->nr_class;
02793 model->label = Malloc(int,n);
02794 for(int i=0;i<n;i++)
02795 fscanf(fp,"%d",&model->label[i]);
02796 }
02797 else if(strcmp(cmd,"probA")==0)
02798 {
02799 int n = model->nr_class * (model->nr_class-1)/2;
02800 model->probA = Malloc(double,n);
02801 for(int i=0;i<n;i++)
02802 fscanf(fp,"%lf",&model->probA[i]);
02803 }
02804 else if(strcmp(cmd,"probB")==0)
02805 {
02806 int n = model->nr_class * (model->nr_class-1)/2;
02807 model->probB = Malloc(double,n);
02808 for(int i=0;i<n;i++)
02809 fscanf(fp,"%lf",&model->probB[i]);
02810 }
02811 else if(strcmp(cmd,"nr_sv")==0)
02812 {
02813 int n = model->nr_class;
02814 model->nSV = Malloc(int,n);
02815 for(int i=0;i<n;i++)
02816 fscanf(fp,"%d",&model->nSV[i]);
02817 }
02818 else if(strcmp(cmd,"SV")==0)
02819 {
02820 while(1)
02821 {
02822 int c = getc(fp);
02823 if(c==EOF || c=='\n') break;
02824 }
02825 break;
02826 }
02827 else
02828 {
02829 fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
02830 free(model->rho);
02831 free(model->label);
02832 free(model->nSV);
02833 free(model);
02834 return NULL;
02835 }
02836 }
02837
02838
02839
02840 int elements = 0;
02841 long pos = ftell(fp);
02842
02843 max_line_len = 1024;
02844 line = Malloc(char,max_line_len);
02845 char *p,*endptr,*idx,*val;
02846
02847 while(readline(fp)!=NULL)
02848 {
02849 p = strtok(line,":");
02850 while(1)
02851 {
02852 p = strtok(NULL,":");
02853 if(p == NULL)
02854 break;
02855 ++elements;
02856 }
02857 }
02858 elements += model->l;
02859
02860 fseek(fp,pos,SEEK_SET);
02861
02862 int m = model->nr_class - 1;
02863 int l = model->l;
02864 model->sv_coef = Malloc(double *,m);
02865 int i;
02866 for(i=0;i<m;i++)
02867 model->sv_coef[i] = Malloc(double,l);
02868 model->SV = Malloc(svm_node*,l);
02869 svm_node *x_space = NULL;
02870 if(l>0) x_space = Malloc(svm_node,elements);
02871
02872 int j=0;
02873 for(i=0;i<l;i++)
02874 {
02875 readline(fp);
02876 model->SV[i] = &x_space[j];
02877
02878 p = strtok(line, " \t");
02879 model->sv_coef[0][i] = strtod(p,&endptr);
02880 for(int k=1;k<m;k++)
02881 {
02882 p = strtok(NULL, " \t");
02883 model->sv_coef[k][i] = strtod(p,&endptr);
02884 }
02885
02886 while(1)
02887 {
02888 idx = strtok(NULL, ":");
02889 val = strtok(NULL, " \t");
02890
02891 if(val == NULL)
02892 break;
02893 x_space[j].index = (int) strtol(idx,&endptr,10);
02894 x_space[j].value = strtod(val,&endptr);
02895
02896 ++j;
02897 }
02898 x_space[j++].index = -1;
02899 }
02900 free(line);
02901
02902 if (ferror(fp) != 0 || fclose(fp) != 0)
02903 return NULL;
02904
02905 model->free_sv = 1;
02906 return model;
02907 }
02908
02909 void svm_free_model_content(svm_model* model_ptr)
02910 {
02911 if(model_ptr->free_sv && model_ptr->l > 0)
02912 free((void *)(model_ptr->SV[0]));
02913 for(int i=0;i<model_ptr->nr_class-1;i++)
02914 free(model_ptr->sv_coef[i]);
02915 free(model_ptr->SV);
02916 free(model_ptr->sv_coef);
02917 free(model_ptr->rho);
02918 free(model_ptr->label);
02919 free(model_ptr->probA);
02920 free(model_ptr->probB);
02921 free(model_ptr->nSV);
02922 free(model_ptr->SVidx);
02923 }
02924
02925 void svm_free_and_destroy_model(svm_model** model_ptr_ptr)
02926 {
02927 svm_model* model_ptr = *model_ptr_ptr;
02928 if(model_ptr != NULL)
02929 {
02930 svm_free_model_content(model_ptr);
02931 free(model_ptr);
02932 }
02933 }
02934
02935 void svm_destroy_model(svm_model* model_ptr)
02936 {
02937 fprintf(stderr,"warning: svm_destroy_model is deprecated and should not be used. Please use svm_free_and_destroy_model(svm_model **model_ptr_ptr)\n");
02938 svm_free_and_destroy_model(&model_ptr);
02939 }
02940
02941 void svm_destroy_param(svm_parameter* param)
02942 {
02943 free(param->weight_label);
02944 free(param->weight);
02945 }
02946
02947 const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
02948 {
02949
02950
02951 int svm_type = param->svm_type;
02952 if(svm_type != C_SVC &&
02953 svm_type != NU_SVC &&
02954 svm_type != ONE_CLASS &&
02955 svm_type != EPSILON_SVR &&
02956 svm_type != NU_SVR)
02957 return "unknown svm type";
02958
02959
02960
02961 int kernel_type = param->kernel_type;
02962 if(kernel_type != LINEAR &&
02963 kernel_type != POLY &&
02964 kernel_type != RBF &&
02965 kernel_type != SIGMOID &&
02966 kernel_type != PRECOMPUTED)
02967 return "unknown kernel type";
02968
02969 if(param->gamma < 0)
02970 return "gamma < 0";
02971
02972 if(param->degree < 0)
02973 return "degree of polynomial kernel < 0";
02974
02975
02976
02977 if(param->cache_size <= 0)
02978 return "cache_size <= 0";
02979
02980 if(param->eps <= 0)
02981 return "eps <= 0";
02982
02983 if(svm_type == C_SVC ||
02984 svm_type == EPSILON_SVR ||
02985 svm_type == NU_SVR)
02986 if(param->C <= 0)
02987 return "C <= 0";
02988
02989 if(svm_type == NU_SVC ||
02990 svm_type == ONE_CLASS ||
02991 svm_type == NU_SVR)
02992 if(param->nu <= 0 || param->nu > 1)
02993 return "nu <= 0 or nu > 1";
02994
02995 if(svm_type == EPSILON_SVR)
02996 if(param->p < 0)
02997 return "p < 0";
02998
02999 if(param->shrinking != 0 &&
03000 param->shrinking != 1)
03001 return "shrinking != 0 and shrinking != 1";
03002
03003 if(param->probability != 0 &&
03004 param->probability != 1)
03005 return "probability != 0 and probability != 1";
03006
03007 if(param->probability == 1 &&
03008 svm_type == ONE_CLASS)
03009 return "one-class SVM probability output not supported yet";
03010
03011
03012
03013
03014 if(svm_type == NU_SVC)
03015 {
03016 int l = prob->l;
03017 int max_nr_class = 16;
03018 int nr_class = 0;
03019 int *label = Malloc(int,max_nr_class);
03020 int *count = Malloc(int,max_nr_class);
03021
03022 int i;
03023 for(i=0;i<l;i++)
03024 {
03025 int this_label = (int)prob->y[i];
03026 int j;
03027 for(j=0;j<nr_class;j++)
03028 if(this_label == label[j])
03029 {
03030 ++count[j];
03031 break;
03032 }
03033 if(j == nr_class)
03034 {
03035 if(nr_class == max_nr_class)
03036 {
03037 max_nr_class *= 2;
03038 label = (int *)realloc(label,max_nr_class*sizeof(int));
03039 count = (int *)realloc(count,max_nr_class*sizeof(int));
03040 }
03041 label[nr_class] = this_label;
03042 count[nr_class] = 1;
03043 ++nr_class;
03044 }
03045 }
03046
03047 for(i=0;i<nr_class;i++)
03048 {
03049 int n1 = count[i];
03050 for(int j=i+1;j<nr_class;j++)
03051 {
03052 int n2 = count[j];
03053 if(param->nu*(n1+n2)/2 > min(n1,n2))
03054 {
03055 free(label);
03056 free(count);
03057 return "specified nu is infeasible";
03058 }
03059 }
03060 }
03061 free(label);
03062 free(count);
03063 }
03064
03065 return NULL;
03066 }
03067
03068 int svm_check_probability_model(const svm_model *model)
03069 {
03070 return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
03071 model->probA!=NULL && model->probB!=NULL) ||
03072 ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
03073 model->probA!=NULL);
03074 }
03075
03076 void svm_set_print_string_function(void (*print_func)(const char *))
03077 {
03078 if(print_func == NULL)
03079 svm_print_string = &print_string_stdout;
03080 else
03081 svm_print_string = print_func;
03082 }