GestureRecognitionToolkit  Version: 0.2.5
The Gesture Recognition Toolkit (GRT) is a cross-platform, open-source, c++ machine learning library for real-time gesture recognition.
libsvm.cpp
1 
2 #include "libsvm.h"
3 
4 namespace LIBSVM {
5 
6 int libsvm_version = LIBSVM_VERSION;
7 typedef float Qfloat;
8 typedef signed char schar;
9 #ifndef min
10 template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
11 #endif
12 #ifndef max
13 template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
14 #endif
15 template <class T> static inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
16 template <class S, class T> static inline void clone(T*& dst, S* src, int n)
17 {
18  dst = new T[n];
19  memcpy((void *)dst,(void *)src,sizeof(T)*n);
20 }
21 static inline double powi(double base, int times)
22 {
23  double tmp = base, ret = 1.0;
24 
25  for(int t=times; t>0; t/=2)
26  {
27  if(t%2==1) ret*=tmp;
28  tmp = tmp * tmp;
29  }
30  return ret;
31 }
32 #define INF HUGE_VAL
33 #define TAU 1e-12
34 #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
35 
36 static void print_string_stdout(const char *s)
37 {
38  fputs(s,stdout);
39  fflush(stdout);
40 }
41 static void (*svm_print_string) (const char *) = &print_string_stdout;
42 #if 1
43 static void info(const char *fmt,...)
44 {
45  char buf[BUFSIZ];
46  va_list ap;
47  va_start(ap,fmt);
48  vsprintf(buf,fmt,ap);
49  va_end(ap);
50  (*svm_print_string)(buf);
51 }
52 #else
53 static void info(const char *fmt,...) {}
54 #endif
55 
56 //
57 // Kernel Cache
58 //
59 // l is the number of total data items
60 // size is the cache size limit in bytes
61 //
62 class Cache
63 {
64 public:
65  Cache(int l,long int size);
66  ~Cache();
67 
68  // request data [0,len)
69  // return some position p where [p,len) need to be filled
70  // (p >= len if nothing needs to be filled)
71  int get_data(const int index, Qfloat **data, int len);
72  void swap_index(int i, int j);
73 private:
74  int l;
75  long int size;
76  struct head_t
77  {
78  head_t *prev, *next; // a circular list
79  Qfloat *data;
80  int len; // data[0,len) is cached in this entry
81  };
82 
83  head_t *head;
84  head_t lru_head;
85  void lru_delete(head_t *h);
86  void lru_insert(head_t *h);
87 };
88 
89 Cache::Cache(int l_,long int size_):l(l_),size(size_)
90 {
91  head = (head_t *)calloc(l,sizeof(head_t)); // initialized to 0
92  size /= sizeof(Qfloat);
93  size -= l * sizeof(head_t) / sizeof(Qfloat);
94  size = std::max(size, 2 * (long int) l); // cache must be large enough for two columns
95  lru_head.next = lru_head.prev = &lru_head;
96 }
97 
98 Cache::~Cache()
99 {
100  for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
101  free(h->data);
102  free(head);
103 }
104 
105 void Cache::lru_delete(head_t *h)
106 {
107  // delete from current location
108  h->prev->next = h->next;
109  h->next->prev = h->prev;
110 }
111 
112 void Cache::lru_insert(head_t *h)
113 {
114  // insert to last position
115  h->next = &lru_head;
116  h->prev = lru_head.prev;
117  h->prev->next = h;
118  h->next->prev = h;
119 }
120 
121 int Cache::get_data(const int index, Qfloat **data, int len)
122 {
123  head_t *h = &head[index];
124  if(h->len) lru_delete(h);
125  int more = len - h->len;
126 
127  if(more > 0)
128  {
129  // free old space
130  while(size < more)
131  {
132  head_t *old = lru_head.next;
133  lru_delete(old);
134  free(old->data);
135  size += old->len;
136  old->data = 0;
137  old->len = 0;
138  }
139 
140  // allocate new space
141  h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
142  size -= more;
143  std::swap(h->len,len);
144  }
145 
146  lru_insert(h);
147  *data = h->data;
148  return len;
149 }
150 
151 void Cache::swap_index(int i, int j)
152 {
153  if(i==j) return;
154 
155  if(head[i].len) lru_delete(&head[i]);
156  if(head[j].len) lru_delete(&head[j]);
157  std::swap(head[i].data,head[j].data);
158  std::swap(head[i].len,head[j].len);
159  if(head[i].len) lru_insert(&head[i]);
160  if(head[j].len) lru_insert(&head[j]);
161 
162  if(i>j) std::swap(i,j);
163  for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
164  {
165  if(h->len > i)
166  {
167  if(h->len > j)
168  std::swap(h->data[i],h->data[j]);
169  else
170  {
171  // give up
172  lru_delete(h);
173  free(h->data);
174  size += h->len;
175  h->data = 0;
176  h->len = 0;
177  }
178  }
179  }
180 }
181 
182 //
183 // Kernel evaluation
184 //
185 // the static method k_function is for doing single kernel evaluation
186 // the constructor of Kernel prepares to calculate the l*l kernel matrix
187 // the member function get_Q is for getting one column from the Q Matrix
188 //
189 class QMatrix {
190 public:
191  virtual Qfloat *get_Q(int column, int len) const = 0;
192  virtual double *get_QD() const = 0;
193  virtual void swap_index(int i, int j) const = 0;
194  virtual ~QMatrix() {}
195 };
196 
197 class Kernel: public QMatrix {
198 public:
199  Kernel(int l, svm_node * const * x, const svm_parameter& param);
200  virtual ~Kernel();
201 
202  static double k_function(const svm_node *x, const svm_node *y,
203  const svm_parameter& param);
204  virtual Qfloat *get_Q(int column, int len) const = 0;
205  virtual double *get_QD() const = 0;
206  virtual void swap_index(int i, int j) const // no so const...
207  {
208  std::swap(x[i],x[j]);
209  if(x_square) std::swap(x_square[i],x_square[j]);
210  }
211 protected:
212 
213  double (Kernel::*kernel_function)(int i, int j) const;
214 
215 private:
216  const svm_node **x;
217  double *x_square;
218 
219  // svm_parameter
220  const int kernel_type;
221  const int degree;
222  const double gamma;
223  const double coef0;
224 
225  static double dot(const svm_node *px, const svm_node *py);
226  double kernel_linear(int i, int j) const
227  {
228  return dot(x[i],x[j]);
229  }
230  double kernel_poly(int i, int j) const
231  {
232  return powi(gamma*dot(x[i],x[j])+coef0,degree);
233  }
234  double kernel_rbf(int i, int j) const
235  {
236  return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
237  }
238  double kernel_sigmoid(int i, int j) const
239  {
240  return tanh(gamma*dot(x[i],x[j])+coef0);
241  }
242  double kernel_precomputed(int i, int j) const
243  {
244  return x[i][(int)(x[j][0].value)].value;
245  }
246 };
247 
248 Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
249 :kernel_type(param.kernel_type), degree(param.degree),
250  gamma(param.gamma), coef0(param.coef0)
251 {
252  switch(kernel_type)
253  {
254  case LINEAR:
255  kernel_function = &Kernel::kernel_linear;
256  break;
257  case POLY:
258  kernel_function = &Kernel::kernel_poly;
259  break;
260  case RBF:
261  kernel_function = &Kernel::kernel_rbf;
262  break;
263  case SIGMOID:
264  kernel_function = &Kernel::kernel_sigmoid;
265  break;
266  case PRECOMPUTED:
267  kernel_function = &Kernel::kernel_precomputed;
268  break;
269  }
270 
271  clone(x,x_,l);
272 
273  if(kernel_type == RBF)
274  {
275  x_square = new double[l];
276  for(int i=0;i<l;i++)
277  x_square[i] = dot(x[i],x[i]);
278  }
279  else
280  x_square = 0;
281 }
282 
283 Kernel::~Kernel()
284 {
285  delete[] x;
286  delete[] x_square;
287 }
288 
289 double Kernel::dot(const svm_node *px, const svm_node *py)
290 {
291  double sum = 0;
292  while(px->index != -1 && py->index != -1)
293  {
294  if(px->index == py->index)
295  {
296  sum += px->value * py->value;
297  ++px;
298  ++py;
299  }
300  else
301  {
302  if(px->index > py->index)
303  ++py;
304  else
305  ++px;
306  }
307  }
308  return sum;
309 }
310 
311 double Kernel::k_function(const svm_node *x, const svm_node *y,
312  const svm_parameter& param)
313 {
314  switch(param.kernel_type)
315  {
316  case LINEAR:
317  return dot(x,y);
318  case POLY:
319  return powi(param.gamma*dot(x,y)+param.coef0,param.degree);
320  case RBF:
321  {
322  double sum = 0;
323  while(x->index != -1 && y->index !=-1)
324  {
325  if(x->index == y->index)
326  {
327  double d = x->value - y->value;
328  sum += d*d;
329  ++x;
330  ++y;
331  }
332  else
333  {
334  if(x->index > y->index)
335  {
336  sum += y->value * y->value;
337  ++y;
338  }
339  else
340  {
341  sum += x->value * x->value;
342  ++x;
343  }
344  }
345  }
346 
347  while(x->index != -1)
348  {
349  sum += x->value * x->value;
350  ++x;
351  }
352 
353  while(y->index != -1)
354  {
355  sum += y->value * y->value;
356  ++y;
357  }
358 
359  return exp(-param.gamma*sum);
360  }
361  case SIGMOID:
362  return tanh(param.gamma*dot(x,y)+param.coef0);
363  case PRECOMPUTED: //x: test (validation), y: SV
364  return x[(int)(y->value)].value;
365  default:
366  return 0; // Unreachable
367  }
368 }
369 
370 // An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
371 // Solves:
372 //
373 // min 0.5(\alpha^T Q \alpha) + p^T \alpha
374 //
375 // y^T \alpha = \delta
376 // y_i = +1 or -1
377 // 0 <= alpha_i <= Cp for y_i = 1
378 // 0 <= alpha_i <= Cn for y_i = -1
379 //
380 // Given:
381 //
382 // Q, p, y, Cp, Cn, and an initial feasible point \alpha
383 // l is the size of vectors and matrices
384 // eps is the stopping tolerance
385 //
386 // solution will be put in \alpha, objective value will be put in obj
387 //
388 class Solver {
389 public:
390  Solver() {};
391  virtual ~Solver() {};
392 
393  struct SolutionInfo {
394  double obj;
395  double rho;
396  double upper_bound_p;
397  double upper_bound_n;
398  double r; // for Solver_NU
399  };
400 
401  void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
402  double *alpha_, double Cp, double Cn, double eps,
403  SolutionInfo* si, int shrinking);
404 protected:
405  int active_size;
406  schar *y;
407  double *G; // gradient of objective function
408  enum { LOWER_BOUND, UPPER_BOUND, FREE };
409  char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
410  double *alpha;
411  const QMatrix *Q;
412  const double *QD;
413  double eps;
414  double Cp,Cn;
415  double *p;
416  int *active_set;
417  double *G_bar; // gradient, if we treat free variables as 0
418  int l;
419  bool unshrink; // XXX
420 
421  double get_C(int i)
422  {
423  return (y[i] > 0)? Cp : Cn;
424  }
425  void update_alpha_status(int i)
426  {
427  if(alpha[i] >= get_C(i))
428  alpha_status[i] = UPPER_BOUND;
429  else if(alpha[i] <= 0)
430  alpha_status[i] = LOWER_BOUND;
431  else alpha_status[i] = FREE;
432  }
433  bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
434  bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
435  bool is_free(int i) { return alpha_status[i] == FREE; }
436  void swap_index(int i, int j);
437  void reconstruct_gradient();
438  virtual int select_working_set(int &i, int &j);
439  virtual double calculate_rho();
440  virtual void do_shrinking();
441 private:
442  bool be_shrunk(int i, double Gmax1, double Gmax2);
443 };
444 
445 void Solver::swap_index(int i, int j)
446 {
447  Q->swap_index(i,j);
448  std::swap(y[i],y[j]);
449  std::swap(G[i],G[j]);
450  std::swap(alpha_status[i],alpha_status[j]);
451  std::swap(alpha[i],alpha[j]);
452  std::swap(p[i],p[j]);
453  std::swap(active_set[i],active_set[j]);
454  std::swap(G_bar[i],G_bar[j]);
455 }
456 
457 void Solver::reconstruct_gradient()
458 {
459  // reconstruct inactive elements of G from G_bar and free variables
460 
461  if(active_size == l) return;
462 
463  int i,j;
464  int nr_free = 0;
465 
466  for(j=active_size;j<l;j++)
467  G[j] = G_bar[j] + p[j];
468 
469  for(j=0;j<active_size;j++)
470  if(is_free(j))
471  nr_free++;
472 
473  if(2*nr_free < active_size)
474  info("\nWARNING: using -h 0 may be faster\n");
475 
476  if (nr_free*l > 2*active_size*(l-active_size))
477  {
478  for(i=active_size;i<l;i++)
479  {
480  const Qfloat *Q_i = Q->get_Q(i,active_size);
481  for(j=0;j<active_size;j++)
482  if(is_free(j))
483  G[i] += alpha[j] * Q_i[j];
484  }
485  }
486  else
487  {
488  for(i=0;i<active_size;i++)
489  if(is_free(i))
490  {
491  const Qfloat *Q_i = Q->get_Q(i,l);
492  double alpha_i = alpha[i];
493  for(j=active_size;j<l;j++)
494  G[j] += alpha_i * Q_i[j];
495  }
496  }
497 }
498 
499 void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
500  double *alpha_, double Cp, double Cn, double eps,
501  SolutionInfo* si, int shrinking)
502 {
503  this->l = l;
504  this->Q = &Q;
505  QD=Q.get_QD();
506  clone(p, p_,l);
507  clone(y, y_,l);
508  clone(alpha,alpha_,l);
509  this->Cp = Cp;
510  this->Cn = Cn;
511  this->eps = eps;
512  unshrink = false;
513 
514  // initialize alpha_status
515  {
516  alpha_status = new char[l];
517  for(int i=0;i<l;i++)
518  update_alpha_status(i);
519  }
520 
521  // initialize active set (for shrinking)
522  {
523  active_set = new int[l];
524  for(int i=0;i<l;i++)
525  active_set[i] = i;
526  active_size = l;
527  }
528 
529  // initialize gradient
530  {
531  G = new double[l];
532  G_bar = new double[l];
533  int i;
534  for(i=0;i<l;i++)
535  {
536  G[i] = p[i];
537  G_bar[i] = 0;
538  }
539  for(i=0;i<l;i++)
540  if(!is_lower_bound(i))
541  {
542  const Qfloat *Q_i = Q.get_Q(i,l);
543  double alpha_i = alpha[i];
544  int j;
545  for(j=0;j<l;j++)
546  G[j] += alpha_i*Q_i[j];
547  if(is_upper_bound(i))
548  for(j=0;j<l;j++)
549  G_bar[j] += get_C(i) * Q_i[j];
550  }
551  }
552 
553  // optimization step
554 
555  int iter = 0;
556  int max_iter = std::max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);
557  int counter = std::min(l,1000)+1;
558 
559  while(iter < max_iter)
560  {
561  // show progress and do shrinking
562 
563  if(--counter == 0)
564  {
565  counter = std::min(l,1000);
566  if(shrinking) do_shrinking();
567  info(".");
568  }
569 
570  int i,j;
571  if(select_working_set(i,j)!=0)
572  {
573  // reconstruct the whole gradient
574  reconstruct_gradient();
575  // reset active set size and check
576  active_size = l;
577  info("*");
578  if(select_working_set(i,j)!=0)
579  break;
580  else
581  counter = 1; // do shrinking next iteration
582  }
583 
584  ++iter;
585 
586  // update alpha[i] and alpha[j], handle bounds carefully
587 
588  const Qfloat *Q_i = Q.get_Q(i,active_size);
589  const Qfloat *Q_j = Q.get_Q(j,active_size);
590 
591  double C_i = get_C(i);
592  double C_j = get_C(j);
593 
594  double old_alpha_i = alpha[i];
595  double old_alpha_j = alpha[j];
596 
597  if(y[i]!=y[j])
598  {
599  double quad_coef = QD[i]+QD[j]+2*Q_i[j];
600  if (quad_coef <= 0)
601  quad_coef = TAU;
602  double delta = (-G[i]-G[j])/quad_coef;
603  double diff = alpha[i] - alpha[j];
604  alpha[i] += delta;
605  alpha[j] += delta;
606 
607  if(diff > 0)
608  {
609  if(alpha[j] < 0)
610  {
611  alpha[j] = 0;
612  alpha[i] = diff;
613  }
614  }
615  else
616  {
617  if(alpha[i] < 0)
618  {
619  alpha[i] = 0;
620  alpha[j] = -diff;
621  }
622  }
623  if(diff > C_i - C_j)
624  {
625  if(alpha[i] > C_i)
626  {
627  alpha[i] = C_i;
628  alpha[j] = C_i - diff;
629  }
630  }
631  else
632  {
633  if(alpha[j] > C_j)
634  {
635  alpha[j] = C_j;
636  alpha[i] = C_j + diff;
637  }
638  }
639  }
640  else
641  {
642  double quad_coef = QD[i]+QD[j]-2*Q_i[j];
643  if (quad_coef <= 0)
644  quad_coef = TAU;
645  double delta = (G[i]-G[j])/quad_coef;
646  double sum = alpha[i] + alpha[j];
647  alpha[i] -= delta;
648  alpha[j] += delta;
649 
650  if(sum > C_i)
651  {
652  if(alpha[i] > C_i)
653  {
654  alpha[i] = C_i;
655  alpha[j] = sum - C_i;
656  }
657  }
658  else
659  {
660  if(alpha[j] < 0)
661  {
662  alpha[j] = 0;
663  alpha[i] = sum;
664  }
665  }
666  if(sum > C_j)
667  {
668  if(alpha[j] > C_j)
669  {
670  alpha[j] = C_j;
671  alpha[i] = sum - C_j;
672  }
673  }
674  else
675  {
676  if(alpha[i] < 0)
677  {
678  alpha[i] = 0;
679  alpha[j] = sum;
680  }
681  }
682  }
683 
684  // update G
685 
686  double delta_alpha_i = alpha[i] - old_alpha_i;
687  double delta_alpha_j = alpha[j] - old_alpha_j;
688 
689  for(int k=0;k<active_size;k++)
690  {
691  G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
692  }
693 
694  // update alpha_status and G_bar
695 
696  {
697  bool ui = is_upper_bound(i);
698  bool uj = is_upper_bound(j);
699  update_alpha_status(i);
700  update_alpha_status(j);
701  int k;
702  if(ui != is_upper_bound(i))
703  {
704  Q_i = Q.get_Q(i,l);
705  if(ui)
706  for(k=0;k<l;k++)
707  G_bar[k] -= C_i * Q_i[k];
708  else
709  for(k=0;k<l;k++)
710  G_bar[k] += C_i * Q_i[k];
711  }
712 
713  if(uj != is_upper_bound(j))
714  {
715  Q_j = Q.get_Q(j,l);
716  if(uj)
717  for(k=0;k<l;k++)
718  G_bar[k] -= C_j * Q_j[k];
719  else
720  for(k=0;k<l;k++)
721  G_bar[k] += C_j * Q_j[k];
722  }
723  }
724  }
725 
726  if(iter >= max_iter)
727  {
728  if(active_size < l)
729  {
730  // reconstruct the whole gradient to calculate objective value
731  reconstruct_gradient();
732  active_size = l;
733  info("*");
734  }
735  info("\nWARNING: reaching max number of iterations");
736  }
737 
738  // calculate rho
739 
740  si->rho = calculate_rho();
741 
742  // calculate objective value
743  {
744  double v = 0;
745  int i;
746  for(i=0;i<l;i++)
747  v += alpha[i] * (G[i] + p[i]);
748 
749  si->obj = v/2;
750  }
751 
752  // put back the solution
753  {
754  for(int i=0;i<l;i++)
755  alpha_[active_set[i]] = alpha[i];
756  }
757 
758  // juggle everything back
759  /*{
760  for(int i=0;i<l;i++)
761  while(active_set[i] != i)
762  swap_index(i,active_set[i]);
763  // or Q.swap_index(i,active_set[i]);
764  }*/
765 
766  si->upper_bound_p = Cp;
767  si->upper_bound_n = Cn;
768 
769  info("\noptimization finished, #iter = %d\n",iter);
770 
771  delete[] p;
772  delete[] y;
773  delete[] alpha;
774  delete[] alpha_status;
775  delete[] active_set;
776  delete[] G;
777  delete[] G_bar;
778 }
779 
780 // return 1 if already optimal, return 0 otherwise
781 int Solver::select_working_set(int &out_i, int &out_j)
782 {
783  // return i,j such that
784  // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
785  // j: minimizes the decrease of obj value
786  // (if quadratic coefficeint <= 0, replace it with tau)
787  // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
788 
789  double Gmax = -INF;
790  double Gmax2 = -INF;
791  int Gmax_idx = -1;
792  int Gmin_idx = -1;
793  double obj_diff_min = INF;
794 
795  for(int t=0;t<active_size;t++)
796  if(y[t]==+1)
797  {
798  if(!is_upper_bound(t))
799  if(-G[t] >= Gmax)
800  {
801  Gmax = -G[t];
802  Gmax_idx = t;
803  }
804  }
805  else
806  {
807  if(!is_lower_bound(t))
808  if(G[t] >= Gmax)
809  {
810  Gmax = G[t];
811  Gmax_idx = t;
812  }
813  }
814 
815  int i = Gmax_idx;
816  const Qfloat *Q_i = NULL;
817  if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
818  Q_i = Q->get_Q(i,active_size);
819 
820  for(int j=0;j<active_size;j++)
821  {
822  if(y[j]==+1)
823  {
824  if (!is_lower_bound(j))
825  {
826  double grad_diff=Gmax+G[j];
827  if (G[j] >= Gmax2)
828  Gmax2 = G[j];
829  if (grad_diff > 0)
830  {
831  double obj_diff;
832  double quad_coef = QD[i]+QD[j]-2.0*y[i]*Q_i[j];
833  if (quad_coef > 0)
834  obj_diff = -(grad_diff*grad_diff)/quad_coef;
835  else
836  obj_diff = -(grad_diff*grad_diff)/TAU;
837 
838  if (obj_diff <= obj_diff_min)
839  {
840  Gmin_idx=j;
841  obj_diff_min = obj_diff;
842  }
843  }
844  }
845  }
846  else
847  {
848  if (!is_upper_bound(j))
849  {
850  double grad_diff= Gmax-G[j];
851  if (-G[j] >= Gmax2)
852  Gmax2 = -G[j];
853  if (grad_diff > 0)
854  {
855  double obj_diff;
856  double quad_coef = QD[i]+QD[j]+2.0*y[i]*Q_i[j];
857  if (quad_coef > 0)
858  obj_diff = -(grad_diff*grad_diff)/quad_coef;
859  else
860  obj_diff = -(grad_diff*grad_diff)/TAU;
861 
862  if (obj_diff <= obj_diff_min)
863  {
864  Gmin_idx=j;
865  obj_diff_min = obj_diff;
866  }
867  }
868  }
869  }
870  }
871 
872  if(Gmax+Gmax2 < eps)
873  return 1;
874 
875  out_i = Gmax_idx;
876  out_j = Gmin_idx;
877  return 0;
878 }
879 
880 bool Solver::be_shrunk(int i, double Gmax1, double Gmax2)
881 {
882  if(is_upper_bound(i))
883  {
884  if(y[i]==+1)
885  return(-G[i] > Gmax1);
886  else
887  return(-G[i] > Gmax2);
888  }
889  else if(is_lower_bound(i))
890  {
891  if(y[i]==+1)
892  return(G[i] > Gmax2);
893  else
894  return(G[i] > Gmax1);
895  }
896  else
897  return(false);
898 }
899 
900 void Solver::do_shrinking()
901 {
902  int i;
903  double Gmax1 = -INF; // max { -y_i * grad(f)_i | i in I_up(\alpha) }
904  double Gmax2 = -INF; // max { y_i * grad(f)_i | i in I_low(\alpha) }
905 
906  // find maximal violating pair first
907  for(i=0;i<active_size;i++)
908  {
909  if(y[i]==+1)
910  {
911  if(!is_upper_bound(i))
912  {
913  if(-G[i] >= Gmax1)
914  Gmax1 = -G[i];
915  }
916  if(!is_lower_bound(i))
917  {
918  if(G[i] >= Gmax2)
919  Gmax2 = G[i];
920  }
921  }
922  else
923  {
924  if(!is_upper_bound(i))
925  {
926  if(-G[i] >= Gmax2)
927  Gmax2 = -G[i];
928  }
929  if(!is_lower_bound(i))
930  {
931  if(G[i] >= Gmax1)
932  Gmax1 = G[i];
933  }
934  }
935  }
936 
937  if(unshrink == false && Gmax1 + Gmax2 <= eps*10)
938  {
939  unshrink = true;
940  reconstruct_gradient();
941  active_size = l;
942  info("*");
943  }
944 
945  for(i=0;i<active_size;i++)
946  if (be_shrunk(i, Gmax1, Gmax2))
947  {
948  active_size--;
949  while (active_size > i)
950  {
951  if (!be_shrunk(active_size, Gmax1, Gmax2))
952  {
953  swap_index(i,active_size);
954  break;
955  }
956  active_size--;
957  }
958  }
959 }
960 
961 double Solver::calculate_rho()
962 {
963  double r;
964  int nr_free = 0;
965  double ub = INF, lb = -INF, sum_free = 0;
966  for(int i=0;i<active_size;i++)
967  {
968  double yG = y[i]*G[i];
969 
970  if(is_upper_bound(i))
971  {
972  if(y[i]==-1)
973  ub = std::min(ub,yG);
974  else
975  lb = std::max(lb,yG);
976  }
977  else if(is_lower_bound(i))
978  {
979  if(y[i]==+1)
980  ub = std::min(ub,yG);
981  else
982  lb = std::max(lb,yG);
983  }
984  else
985  {
986  ++nr_free;
987  sum_free += yG;
988  }
989  }
990 
991  if(nr_free>0)
992  r = sum_free/nr_free;
993  else
994  r = (ub+lb)/2;
995 
996  return r;
997 }
998 
999 //
1000 // Solver for nu-svm classification and regression
1001 //
1002 // additional constraint: e^T \alpha = constant
1003 //
1004 class Solver_NU : public Solver
1005 {
1006 public:
1007  Solver_NU() {}
1008  void Solve(int l, const QMatrix& Q, const double *p, const schar *y,
1009  double *alpha, double Cp, double Cn, double eps,
1010  SolutionInfo* si, int shrinking)
1011  {
1012  this->si = si;
1013  Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking);
1014  }
1015 private:
1016  SolutionInfo *si;
1017  int select_working_set(int &i, int &j);
1018  double calculate_rho();
1019  bool be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4);
1020  void do_shrinking();
1021 };
1022 
1023 // return 1 if already optimal, return 0 otherwise
1024 int Solver_NU::select_working_set(int &out_i, int &out_j)
1025 {
1026  // return i,j such that y_i = y_j and
1027  // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
1028  // j: minimizes the decrease of obj value
1029  // (if quadratic coefficeint <= 0, replace it with tau)
1030  // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
1031 
1032  double Gmaxp = -INF;
1033  double Gmaxp2 = -INF;
1034  int Gmaxp_idx = -1;
1035 
1036  double Gmaxn = -INF;
1037  double Gmaxn2 = -INF;
1038  int Gmaxn_idx = -1;
1039 
1040  int Gmin_idx = -1;
1041  double obj_diff_min = INF;
1042 
1043  for(int t=0;t<active_size;t++)
1044  if(y[t]==+1)
1045  {
1046  if(!is_upper_bound(t))
1047  if(-G[t] >= Gmaxp)
1048  {
1049  Gmaxp = -G[t];
1050  Gmaxp_idx = t;
1051  }
1052  }
1053  else
1054  {
1055  if(!is_lower_bound(t))
1056  if(G[t] >= Gmaxn)
1057  {
1058  Gmaxn = G[t];
1059  Gmaxn_idx = t;
1060  }
1061  }
1062 
1063  int ip = Gmaxp_idx;
1064  int in = Gmaxn_idx;
1065  const Qfloat *Q_ip = NULL;
1066  const Qfloat *Q_in = NULL;
1067  if(ip != -1) // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1
1068  Q_ip = Q->get_Q(ip,active_size);
1069  if(in != -1)
1070  Q_in = Q->get_Q(in,active_size);
1071 
1072  for(int j=0;j<active_size;j++)
1073  {
1074  if(y[j]==+1)
1075  {
1076  if (!is_lower_bound(j))
1077  {
1078  double grad_diff=Gmaxp+G[j];
1079  if (G[j] >= Gmaxp2)
1080  Gmaxp2 = G[j];
1081  if (grad_diff > 0)
1082  {
1083  double obj_diff;
1084  double quad_coef = QD[ip]+QD[j]-2*Q_ip[j];
1085  if (quad_coef > 0)
1086  obj_diff = -(grad_diff*grad_diff)/quad_coef;
1087  else
1088  obj_diff = -(grad_diff*grad_diff)/TAU;
1089 
1090  if (obj_diff <= obj_diff_min)
1091  {
1092  Gmin_idx=j;
1093  obj_diff_min = obj_diff;
1094  }
1095  }
1096  }
1097  }
1098  else
1099  {
1100  if (!is_upper_bound(j))
1101  {
1102  double grad_diff=Gmaxn-G[j];
1103  if (-G[j] >= Gmaxn2)
1104  Gmaxn2 = -G[j];
1105  if (grad_diff > 0)
1106  {
1107  double obj_diff;
1108  double quad_coef = QD[in]+QD[j]-2*Q_in[j];
1109  if (quad_coef > 0)
1110  obj_diff = -(grad_diff*grad_diff)/quad_coef;
1111  else
1112  obj_diff = -(grad_diff*grad_diff)/TAU;
1113 
1114  if (obj_diff <= obj_diff_min)
1115  {
1116  Gmin_idx=j;
1117  obj_diff_min = obj_diff;
1118  }
1119  }
1120  }
1121  }
1122  }
1123 
1124  if(std::max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
1125  return 1;
1126 
1127  if (y[Gmin_idx] == +1)
1128  out_i = Gmaxp_idx;
1129  else
1130  out_i = Gmaxn_idx;
1131  out_j = Gmin_idx;
1132 
1133  return 0;
1134 }
1135 
1136 bool Solver_NU::be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)
1137 {
1138  if(is_upper_bound(i))
1139  {
1140  if(y[i]==+1)
1141  return(-G[i] > Gmax1);
1142  else
1143  return(-G[i] > Gmax4);
1144  }
1145  else if(is_lower_bound(i))
1146  {
1147  if(y[i]==+1)
1148  return(G[i] > Gmax2);
1149  else
1150  return(G[i] > Gmax3);
1151  }
1152  else
1153  return(false);
1154 }
1155 
1156 void Solver_NU::do_shrinking()
1157 {
1158  double Gmax1 = -INF; // max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }
1159  double Gmax2 = -INF; // max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }
1160  double Gmax3 = -INF; // max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
1161  double Gmax4 = -INF; // max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }
1162 
1163  // find maximal violating pair first
1164  int i;
1165  for(i=0;i<active_size;i++)
1166  {
1167  if(!is_upper_bound(i))
1168  {
1169  if(y[i]==+1)
1170  {
1171  if(-G[i] > Gmax1) Gmax1 = -G[i];
1172  }
1173  else if(-G[i] > Gmax4) Gmax4 = -G[i];
1174  }
1175  if(!is_lower_bound(i))
1176  {
1177  if(y[i]==+1)
1178  {
1179  if(G[i] > Gmax2) Gmax2 = G[i];
1180  }
1181  else if(G[i] > Gmax3) Gmax3 = G[i];
1182  }
1183  }
1184 
1185  if(unshrink == false && std::max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10)
1186  {
1187  unshrink = true;
1188  reconstruct_gradient();
1189  active_size = l;
1190  }
1191 
1192  for(i=0;i<active_size;i++)
1193  if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
1194  {
1195  active_size--;
1196  while (active_size > i)
1197  {
1198  if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
1199  {
1200  swap_index(i,active_size);
1201  break;
1202  }
1203  active_size--;
1204  }
1205  }
1206 }
1207 
1208 double Solver_NU::calculate_rho()
1209 {
1210  int nr_free1 = 0,nr_free2 = 0;
1211  double ub1 = INF, ub2 = INF;
1212  double lb1 = -INF, lb2 = -INF;
1213  double sum_free1 = 0, sum_free2 = 0;
1214 
1215  for(int i=0;i<active_size;i++)
1216  {
1217  if(y[i]==+1)
1218  {
1219  if(is_upper_bound(i))
1220  lb1 = std::max(lb1,G[i]);
1221  else if(is_lower_bound(i))
1222  ub1 = std::min(ub1,G[i]);
1223  else
1224  {
1225  ++nr_free1;
1226  sum_free1 += G[i];
1227  }
1228  }
1229  else
1230  {
1231  if(is_upper_bound(i))
1232  lb2 = std::max(lb2,G[i]);
1233  else if(is_lower_bound(i))
1234  ub2 = std::min(ub2,G[i]);
1235  else
1236  {
1237  ++nr_free2;
1238  sum_free2 += G[i];
1239  }
1240  }
1241  }
1242 
1243  double r1,r2;
1244  if(nr_free1 > 0)
1245  r1 = sum_free1/nr_free1;
1246  else
1247  r1 = (ub1+lb1)/2;
1248 
1249  if(nr_free2 > 0)
1250  r2 = sum_free2/nr_free2;
1251  else
1252  r2 = (ub2+lb2)/2;
1253 
1254  si->r = (r1+r2)/2;
1255  return (r1-r2)/2;
1256 }
1257 
1258 //
1259 // Q matrices for various formulations
1260 //
1261 class SVC_Q: public Kernel
1262 {
1263 public:
1264  SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
1265  :Kernel(prob.l, prob.x, param)
1266  {
1267  clone(y,y_,prob.l);
1268  cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
1269  QD = new double[prob.l];
1270  for(int i=0;i<prob.l;i++)
1271  QD[i] = (this->*kernel_function)(i,i);
1272  }
1273 
1274  Qfloat *get_Q(int i, int len) const
1275  {
1276  Qfloat *data;
1277  int start, j;
1278  if((start = cache->get_data(i,&data,len)) < len)
1279  {
1280  for(j=start;j<len;j++)
1281  data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
1282  }
1283  return data;
1284  }
1285 
1286  double *get_QD() const
1287  {
1288  return QD;
1289  }
1290 
1291  void swap_index(int i, int j) const
1292  {
1293  cache->swap_index(i,j);
1294  Kernel::swap_index(i,j);
1295  std::swap(y[i],y[j]);
1296  std::swap(QD[i],QD[j]);
1297  }
1298 
1299  ~SVC_Q()
1300  {
1301  delete[] y;
1302  delete cache;
1303  delete[] QD;
1304  }
1305 private:
1306  schar *y;
1307  Cache *cache;
1308  double *QD;
1309 };
1310 
1311 class ONE_CLASS_Q: public Kernel
1312 {
1313 public:
1314  ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
1315  :Kernel(prob.l, prob.x, param)
1316  {
1317  cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
1318  QD = new double[prob.l];
1319  for(int i=0;i<prob.l;i++)
1320  QD[i] = (this->*kernel_function)(i,i);
1321  }
1322 
1323  Qfloat *get_Q(int i, int len) const
1324  {
1325  Qfloat *data;
1326  int start, j;
1327  if((start = cache->get_data(i,&data,len)) < len)
1328  {
1329  for(j=start;j<len;j++)
1330  data[j] = (Qfloat)(this->*kernel_function)(i,j);
1331  }
1332  return data;
1333  }
1334 
1335  double *get_QD() const
1336  {
1337  return QD;
1338  }
1339 
1340  void swap_index(int i, int j) const
1341  {
1342  cache->swap_index(i,j);
1343  Kernel::swap_index(i,j);
1344  std::swap(QD[i],QD[j]);
1345  }
1346 
1347  ~ONE_CLASS_Q()
1348  {
1349  delete cache;
1350  delete[] QD;
1351  }
1352 private:
1353  Cache *cache;
1354  double *QD;
1355 };
1356 
1357 class SVR_Q: public Kernel
1358 {
1359 public:
1360  SVR_Q(const svm_problem& prob, const svm_parameter& param)
1361  :Kernel(prob.l, prob.x, param)
1362  {
1363  l = prob.l;
1364  cache = new Cache(l,(long int)(param.cache_size*(1<<20)));
1365  QD = new double[2*l];
1366  sign = new schar[2*l];
1367  index = new int[2*l];
1368  for(int k=0;k<l;k++)
1369  {
1370  sign[k] = 1;
1371  sign[k+l] = -1;
1372  index[k] = k;
1373  index[k+l] = k;
1374  QD[k] = (this->*kernel_function)(k,k);
1375  QD[k+l] = QD[k];
1376  }
1377  buffer[0] = new Qfloat[2*l];
1378  buffer[1] = new Qfloat[2*l];
1379  next_buffer = 0;
1380  }
1381 
1382  void swap_index(int i, int j) const
1383  {
1384  std::swap(sign[i],sign[j]);
1385  std::swap(index[i],index[j]);
1386  std::swap(QD[i],QD[j]);
1387  }
1388 
1389  Qfloat *get_Q(int i, int len) const
1390  {
1391  Qfloat *data;
1392  int j, real_i = index[i];
1393  if(cache->get_data(real_i,&data,l) < l)
1394  {
1395  for(j=0;j<l;j++)
1396  data[j] = (Qfloat)(this->*kernel_function)(real_i,j);
1397  }
1398 
1399  // reorder and copy
1400  Qfloat *buf = buffer[next_buffer];
1401  next_buffer = 1 - next_buffer;
1402  schar si = sign[i];
1403  for(j=0;j<len;j++)
1404  buf[j] = (Qfloat) si * (Qfloat) sign[j] * data[index[j]];
1405  return buf;
1406  }
1407 
1408  double *get_QD() const
1409  {
1410  return QD;
1411  }
1412 
1413  ~SVR_Q()
1414  {
1415  delete cache;
1416  delete[] sign;
1417  delete[] index;
1418  delete[] buffer[0];
1419  delete[] buffer[1];
1420  delete[] QD;
1421  }
1422 private:
1423  int l;
1424  Cache *cache;
1425  schar *sign;
1426  int *index;
1427  mutable int next_buffer;
1428  Qfloat *buffer[2];
1429  double *QD;
1430 };
1431 
1432 //
1433 // construct and solve various formulations
1434 //
1435 static void solve_c_svc(
1436  const svm_problem *prob, const svm_parameter* param,
1437  double *alpha, Solver::SolutionInfo* si, double Cp, double Cn)
1438 {
1439  int l = prob->l;
1440  double *minus_ones = new double[l];
1441  schar *y = new schar[l];
1442 
1443  int i;
1444 
1445  for(i=0;i<l;i++)
1446  {
1447  alpha[i] = 0;
1448  minus_ones[i] = -1;
1449  if(prob->y[i] > 0) y[i] = +1; else y[i] = -1;
1450  }
1451 
1452  Solver s;
1453  s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
1454  alpha, Cp, Cn, param->eps, si, param->shrinking);
1455 
1456  double sum_alpha=0;
1457  for(i=0;i<l;i++)
1458  sum_alpha += alpha[i];
1459 
1460  if (Cp==Cn)
1461  info("nu = %f\n", sum_alpha/(Cp*prob->l));
1462 
1463  for(i=0;i<l;i++)
1464  alpha[i] *= y[i];
1465 
1466  delete[] minus_ones;
1467  delete[] y;
1468 }
1469 
1470 static void solve_nu_svc(
1471  const svm_problem *prob, const svm_parameter *param,
1472  double *alpha, Solver::SolutionInfo* si)
1473 {
1474  int i;
1475  int l = prob->l;
1476  double nu = param->nu;
1477 
1478  schar *y = new schar[l];
1479 
1480  for(i=0;i<l;i++)
1481  if(prob->y[i]>0)
1482  y[i] = +1;
1483  else
1484  y[i] = -1;
1485 
1486  double sum_pos = nu*l/2;
1487  double sum_neg = nu*l/2;
1488 
1489  for(i=0;i<l;i++)
1490  if(y[i] == +1)
1491  {
1492  alpha[i] = std::min(1.0,sum_pos);
1493  sum_pos -= alpha[i];
1494  }
1495  else
1496  {
1497  alpha[i] = std::min(1.0,sum_neg);
1498  sum_neg -= alpha[i];
1499  }
1500 
1501  double *zeros = new double[l];
1502 
1503  for(i=0;i<l;i++)
1504  zeros[i] = 0;
1505 
1506  Solver_NU s;
1507  s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
1508  alpha, 1.0, 1.0, param->eps, si, param->shrinking);
1509  double r = si->r;
1510 
1511  info("C = %f\n",1/r);
1512 
1513  for(i=0;i<l;i++)
1514  alpha[i] *= y[i]/r;
1515 
1516  si->rho /= r;
1517  si->obj /= (r*r);
1518  si->upper_bound_p = 1/r;
1519  si->upper_bound_n = 1/r;
1520 
1521  delete[] y;
1522  delete[] zeros;
1523 }
1524 
1525 static void solve_one_class(
1526  const svm_problem *prob, const svm_parameter *param,
1527  double *alpha, Solver::SolutionInfo* si)
1528 {
1529  int l = prob->l;
1530  double *zeros = new double[l];
1531  schar *ones = new schar[l];
1532  int i;
1533 
1534  int n = (int)(param->nu*prob->l); // # of alpha's at upper bound
1535 
1536  for(i=0;i<n;i++)
1537  alpha[i] = 1;
1538  if(n<prob->l)
1539  alpha[n] = param->nu * prob->l - n;
1540  for(i=n+1;i<l;i++)
1541  alpha[i] = 0;
1542 
1543  for(i=0;i<l;i++)
1544  {
1545  zeros[i] = 0;
1546  ones[i] = 1;
1547  }
1548 
1549  Solver s;
1550  s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
1551  alpha, 1.0, 1.0, param->eps, si, param->shrinking);
1552 
1553  delete[] zeros;
1554  delete[] ones;
1555 }
1556 
1557 static void solve_epsilon_svr(
1558  const svm_problem *prob, const svm_parameter *param,
1559  double *alpha, Solver::SolutionInfo* si)
1560 {
1561  int l = prob->l;
1562  double *alpha2 = new double[2*l];
1563  double *linear_term = new double[2*l];
1564  schar *y = new schar[2*l];
1565  int i;
1566 
1567  for(i=0;i<l;i++)
1568  {
1569  alpha2[i] = 0;
1570  linear_term[i] = param->p - prob->y[i];
1571  y[i] = 1;
1572 
1573  alpha2[i+l] = 0;
1574  linear_term[i+l] = param->p + prob->y[i];
1575  y[i+l] = -1;
1576  }
1577 
1578  Solver s;
1579  s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
1580  alpha2, param->C, param->C, param->eps, si, param->shrinking);
1581 
1582  double sum_alpha = 0;
1583  for(i=0;i<l;i++)
1584  {
1585  alpha[i] = alpha2[i] - alpha2[i+l];
1586  sum_alpha += fabs(alpha[i]);
1587  }
1588  info("nu = %f\n",sum_alpha/(param->C*l));
1589 
1590  delete[] alpha2;
1591  delete[] linear_term;
1592  delete[] y;
1593 }
1594 
1595 static void solve_nu_svr(
1596  const svm_problem *prob, const svm_parameter *param,
1597  double *alpha, Solver::SolutionInfo* si)
1598 {
1599  int l = prob->l;
1600  double C = param->C;
1601  double *alpha2 = new double[2*l];
1602  double *linear_term = new double[2*l];
1603  schar *y = new schar[2*l];
1604  int i;
1605 
1606  double sum = C * param->nu * l / 2;
1607  for(i=0;i<l;i++)
1608  {
1609  alpha2[i] = alpha2[i+l] = std::min(sum,C);
1610  sum -= alpha2[i];
1611 
1612  linear_term[i] = - prob->y[i];
1613  y[i] = 1;
1614 
1615  linear_term[i+l] = prob->y[i];
1616  y[i+l] = -1;
1617  }
1618 
1619  Solver_NU s;
1620  s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
1621  alpha2, C, C, param->eps, si, param->shrinking);
1622 
1623  info("epsilon = %f\n",-si->r);
1624 
1625  for(i=0;i<l;i++)
1626  alpha[i] = alpha2[i] - alpha2[i+l];
1627 
1628  delete[] alpha2;
1629  delete[] linear_term;
1630  delete[] y;
1631 }
1632 
1633 //
1634 // decision_function
1635 //
1637 {
1638  double *alpha;
1639  double rho;
1640 };
1641 
1642 static decision_function svm_train_one(
1643  const svm_problem *prob, const svm_parameter *param,
1644  double Cp, double Cn)
1645 {
1646  double *alpha = Malloc(double,prob->l);
1648  switch(param->svm_type)
1649  {
1650  case C_SVC:
1651  solve_c_svc(prob,param,alpha,&si,Cp,Cn);
1652  break;
1653  case NU_SVC:
1654  solve_nu_svc(prob,param,alpha,&si);
1655  break;
1656  case ONE_CLASS:
1657  solve_one_class(prob,param,alpha,&si);
1658  break;
1659  case EPSILON_SVR:
1660  solve_epsilon_svr(prob,param,alpha,&si);
1661  break;
1662  case NU_SVR:
1663  solve_nu_svr(prob,param,alpha,&si);
1664  break;
1665  }
1666 
1667  info("obj = %f, rho = %f\n",si.obj,si.rho);
1668 
1669  // output SVs
1670 
1671  int nSV = 0;
1672  int nBSV = 0;
1673  for(int i=0;i<prob->l;i++)
1674  {
1675  if(fabs(alpha[i]) > 0)
1676  {
1677  ++nSV;
1678  if(prob->y[i] > 0)
1679  {
1680  if(fabs(alpha[i]) >= si.upper_bound_p)
1681  ++nBSV;
1682  }
1683  else
1684  {
1685  if(fabs(alpha[i]) >= si.upper_bound_n)
1686  ++nBSV;
1687  }
1688  }
1689  }
1690 
1691  info("nSV = %d, nBSV = %d\n",nSV,nBSV);
1692 
1694  f.alpha = alpha;
1695  f.rho = si.rho;
1696  return f;
1697 }
1698 
1699 // Platt's binary SVM Probablistic Output: an improvement from Lin et al.
1700 static void sigmoid_train(
1701  int l, const double *dec_values, const double *labels,
1702  double& A, double& B)
1703 {
1704  double prior1=0, prior0 = 0;
1705  int i;
1706 
1707  for (i=0;i<l;i++)
1708  if (labels[i] > 0) prior1+=1;
1709  else prior0+=1;
1710 
1711  int max_iter=100; // Maximal number of iterations
1712  double min_step=1e-10; // Minimal step taken in line search
1713  double sigma=1e-12; // For numerically strict PD of Hessian
1714  double eps=1e-5;
1715  double hiTarget=(prior1+1.0)/(prior1+2.0);
1716  double loTarget=1/(prior0+2.0);
1717  double *t=Malloc(double,l);
1718  double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
1719  double newA,newB,newf,d1,d2;
1720  int iter;
1721 
1722  // Initial Point and Initial Fun Value
1723  A=0.0; B=log((prior0+1.0)/(prior1+1.0));
1724  double fval = 0.0;
1725 
1726  for (i=0;i<l;i++)
1727  {
1728  if (labels[i]>0) t[i]=hiTarget;
1729  else t[i]=loTarget;
1730  fApB = dec_values[i]*A+B;
1731  if (fApB>=0)
1732  fval += t[i]*fApB + log(1+exp(-fApB));
1733  else
1734  fval += (t[i] - 1)*fApB +log(1+exp(fApB));
1735  }
1736  for (iter=0;iter<max_iter;iter++)
1737  {
1738  // Update Gradient and Hessian (use H' = H + sigma I)
1739  h11=sigma; // numerically ensures strict PD
1740  h22=sigma;
1741  h21=0.0;g1=0.0;g2=0.0;
1742  for (i=0;i<l;i++)
1743  {
1744  fApB = dec_values[i]*A+B;
1745  if (fApB >= 0)
1746  {
1747  p=exp(-fApB)/(1.0+exp(-fApB));
1748  q=1.0/(1.0+exp(-fApB));
1749  }
1750  else
1751  {
1752  p=1.0/(1.0+exp(fApB));
1753  q=exp(fApB)/(1.0+exp(fApB));
1754  }
1755  d2=p*q;
1756  h11+=dec_values[i]*dec_values[i]*d2;
1757  h22+=d2;
1758  h21+=dec_values[i]*d2;
1759  d1=t[i]-p;
1760  g1+=dec_values[i]*d1;
1761  g2+=d1;
1762  }
1763 
1764  // Stopping Criteria
1765  if (fabs(g1)<eps && fabs(g2)<eps)
1766  break;
1767 
1768  // Finding Newton direction: -inv(H') * g
1769  det=h11*h22-h21*h21;
1770  dA=-(h22*g1 - h21 * g2) / det;
1771  dB=-(-h21*g1+ h11 * g2) / det;
1772  gd=g1*dA+g2*dB;
1773 
1774 
1775  stepsize = 1; // Line Search
1776  while (stepsize >= min_step)
1777  {
1778  newA = A + stepsize * dA;
1779  newB = B + stepsize * dB;
1780 
1781  // New function value
1782  newf = 0.0;
1783  for (i=0;i<l;i++)
1784  {
1785  fApB = dec_values[i]*newA+newB;
1786  if (fApB >= 0)
1787  newf += t[i]*fApB + log(1+exp(-fApB));
1788  else
1789  newf += (t[i] - 1)*fApB +log(1+exp(fApB));
1790  }
1791  // Check sufficient decrease
1792  if (newf<fval+0.0001*stepsize*gd)
1793  {
1794  A=newA;B=newB;fval=newf;
1795  break;
1796  }
1797  else
1798  stepsize = stepsize / 2.0;
1799  }
1800 
1801  if (stepsize < min_step)
1802  {
1803  info("Line search fails in two-class probability estimates\n");
1804  break;
1805  }
1806  }
1807 
1808  if (iter>=max_iter)
1809  info("Reaching maximal iterations in two-class probability estimates\n");
1810  free(t);
1811 }
1812 
1813 static double sigmoid_predict(double decision_value, double A, double B)
1814 {
1815  double fApB = decision_value*A+B;
1816  // 1-p used later; avoid catastrophic cancellation
1817  if (fApB >= 0)
1818  return exp(-fApB)/(1.0+exp(-fApB));
1819  else
1820  return 1.0/(1+exp(fApB)) ;
1821 }
1822 
1823 // Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
1824 static void multiclass_probability(int k, double **r, double *p)
1825 {
1826  int t,j;
1827  int iter = 0, max_iter=std::max(100,k);
1828  double **Q=Malloc(double *,k);
1829  double *Qp=Malloc(double,k);
1830  double pQp, eps=0.005/k;
1831 
1832  for (t=0;t<k;t++)
1833  {
1834  p[t]=1.0/k; // Valid if k = 1
1835  Q[t]=Malloc(double,k);
1836  Q[t][t]=0;
1837  for (j=0;j<t;j++)
1838  {
1839  Q[t][t]+=r[j][t]*r[j][t];
1840  Q[t][j]=Q[j][t];
1841  }
1842  for (j=t+1;j<k;j++)
1843  {
1844  Q[t][t]+=r[j][t]*r[j][t];
1845  Q[t][j]=-r[j][t]*r[t][j];
1846  }
1847  }
1848  for (iter=0;iter<max_iter;iter++)
1849  {
1850  // stopping condition, recalculate QP,pQP for numerical accuracy
1851  pQp=0;
1852  for (t=0;t<k;t++)
1853  {
1854  Qp[t]=0;
1855  for (j=0;j<k;j++)
1856  Qp[t]+=Q[t][j]*p[j];
1857  pQp+=p[t]*Qp[t];
1858  }
1859  double max_error=0;
1860  for (t=0;t<k;t++)
1861  {
1862  double error=fabs(Qp[t]-pQp);
1863  if (error>max_error)
1864  max_error=error;
1865  }
1866  if (max_error<eps) break;
1867 
1868  for (t=0;t<k;t++)
1869  {
1870  double diff=(-Qp[t]+pQp)/Q[t][t];
1871  p[t]+=diff;
1872  pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
1873  for (j=0;j<k;j++)
1874  {
1875  Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
1876  p[j]/=(1+diff);
1877  }
1878  }
1879  }
1880  if (iter>=max_iter)
1881  info("Exceeds max_iter in multiclass_prob\n");
1882  for(t=0;t<k;t++) free(Q[t]);
1883  free(Q);
1884  free(Qp);
1885 }
1886 
1887 // Cross-validation decision values for probability estimates
1888 static void svm_binary_svc_probability(
1889  const svm_problem *prob, const svm_parameter *param,
1890  double Cp, double Cn, double& probA, double& probB)
1891 {
1892  int i;
1893  int nr_fold = 5;
1894  int *perm = Malloc(int,prob->l);
1895  double *dec_values = Malloc(double,prob->l);
1896 
1897  // random shuffle
1898  for(i=0;i<prob->l;i++) perm[i]=i;
1899  for(i=0;i<prob->l;i++)
1900  {
1901  int j = i+rand()%(prob->l-i);
1902  std::swap(perm[i],perm[j]);
1903  }
1904  for(i=0;i<nr_fold;i++)
1905  {
1906  int begin = i*prob->l/nr_fold;
1907  int end = (i+1)*prob->l/nr_fold;
1908  int j,k;
1909  struct svm_problem subprob;
1910 
1911  subprob.l = prob->l-(end-begin);
1912  subprob.x = Malloc(struct svm_node*,subprob.l);
1913  subprob.y = Malloc(double,subprob.l);
1914 
1915  k=0;
1916  for(j=0;j<begin;j++)
1917  {
1918  subprob.x[k] = prob->x[perm[j]];
1919  subprob.y[k] = prob->y[perm[j]];
1920  ++k;
1921  }
1922  for(j=end;j<prob->l;j++)
1923  {
1924  subprob.x[k] = prob->x[perm[j]];
1925  subprob.y[k] = prob->y[perm[j]];
1926  ++k;
1927  }
1928  int p_count=0,n_count=0;
1929  for(j=0;j<k;j++)
1930  if(subprob.y[j]>0)
1931  p_count++;
1932  else
1933  n_count++;
1934 
1935  if(p_count==0 && n_count==0)
1936  for(j=begin;j<end;j++)
1937  dec_values[perm[j]] = 0;
1938  else if(p_count > 0 && n_count == 0)
1939  for(j=begin;j<end;j++)
1940  dec_values[perm[j]] = 1;
1941  else if(p_count == 0 && n_count > 0)
1942  for(j=begin;j<end;j++)
1943  dec_values[perm[j]] = -1;
1944  else
1945  {
1946  svm_parameter subparam = *param;
1947  subparam.probability=0;
1948  subparam.C=1.0;
1949  subparam.nr_weight=2;
1950  subparam.weight_label = Malloc(int,2);
1951  subparam.weight = Malloc(double,2);
1952  subparam.weight_label[0]=+1;
1953  subparam.weight_label[1]=-1;
1954  subparam.weight[0]=Cp;
1955  subparam.weight[1]=Cn;
1956  struct svm_model *submodel = svm_train(&subprob,&subparam);
1957  for(j=begin;j<end;j++)
1958  {
1959  svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]]));
1960  // ensure +1 -1 order; reason not using CV subroutine
1961  dec_values[perm[j]] *= submodel->label[0];
1962  }
1963  svm_free_and_destroy_model(&submodel);
1964  svm_destroy_param(&subparam);
1965  }
1966  free(subprob.x);
1967  free(subprob.y);
1968  }
1969  sigmoid_train(prob->l,dec_values,prob->y,probA,probB);
1970  free(dec_values);
1971  free(perm);
1972 }
1973 
1974 // Return parameter of a Laplace distribution
1975 static double svm_svr_probability(
1976  const svm_problem *prob, const svm_parameter *param)
1977 {
1978  int i;
1979  int nr_fold = 5;
1980  double *ymv = Malloc(double,prob->l);
1981  double mae = 0;
1982 
1983  svm_parameter newparam = *param;
1984  newparam.probability = 0;
1985  svm_cross_validation(prob,&newparam,nr_fold,ymv);
1986  for(i=0;i<prob->l;i++)
1987  {
1988  ymv[i]=prob->y[i]-ymv[i];
1989  mae += fabs(ymv[i]);
1990  }
1991  mae /= prob->l;
1992  double std=sqrt(2*mae*mae);
1993  int count=0;
1994  mae=0;
1995  for(i=0;i<prob->l;i++)
1996  if (fabs(ymv[i]) > 5*std)
1997  count=count+1;
1998  else
1999  mae+=fabs(ymv[i]);
2000  mae /= (prob->l-count);
2001  info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
2002  free(ymv);
2003  return mae;
2004 }
2005 
2006 
2007 // label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
2008 // perm, length l, must be allocated before calling this subroutine
2009 static void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)
2010 {
2011  int l = prob->l;
2012  int max_nr_class = 16;
2013  int nr_class = 0;
2014  int *label = Malloc(int,max_nr_class);
2015  int *count = Malloc(int,max_nr_class);
2016  int *data_label = Malloc(int,l);
2017  int i;
2018 
2019  for(i=0;i<l;i++)
2020  {
2021  int this_label = (int)prob->y[i];
2022  int j;
2023  for(j=0;j<nr_class;j++)
2024  {
2025  if(this_label == label[j])
2026  {
2027  ++count[j];
2028  break;
2029  }
2030  }
2031  data_label[i] = j;
2032  if(j == nr_class)
2033  {
2034  if(nr_class == max_nr_class)
2035  {
2036  max_nr_class *= 2;
2037  label = (int *)realloc(label,max_nr_class*sizeof(int));
2038  count = (int *)realloc(count,max_nr_class*sizeof(int));
2039  }
2040  label[nr_class] = this_label;
2041  count[nr_class] = 1;
2042  ++nr_class;
2043  }
2044  }
2045 
2046  int *start = Malloc(int,nr_class);
2047  start[0] = 0;
2048  for(i=1;i<nr_class;i++)
2049  start[i] = start[i-1]+count[i-1];
2050  for(i=0;i<l;i++)
2051  {
2052  perm[start[data_label[i]]] = i;
2053  ++start[data_label[i]];
2054  }
2055  start[0] = 0;
2056  for(i=1;i<nr_class;i++)
2057  start[i] = start[i-1]+count[i-1];
2058 
2059  *nr_class_ret = nr_class;
2060  *label_ret = label;
2061  *start_ret = start;
2062  *count_ret = count;
2063  free(data_label);
2064 }
2065 
2066 //
2067 // Interface functions
2068 //
2069 svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
2070 {
2071  svm_model *model = Malloc(svm_model,1);
2072  model->param = *param;
2073  model->free_sv = 0; // XXX
2074 
2075  if(param->svm_type == ONE_CLASS ||
2076  param->svm_type == EPSILON_SVR ||
2077  param->svm_type == NU_SVR)
2078  {
2079  // regression or one-class-svm
2080  model->nr_class = 2;
2081  model->label = NULL;
2082  model->nSV = NULL;
2083  model->probA = NULL; model->probB = NULL;
2084  model->sv_coef = Malloc(double *,1);
2085 
2086  if(param->probability &&
2087  (param->svm_type == EPSILON_SVR ||
2088  param->svm_type == NU_SVR))
2089  {
2090  model->probA = Malloc(double,1);
2091  model->probA[0] = svm_svr_probability(prob,param);
2092  }
2093 
2094  decision_function f = svm_train_one(prob,param,0,0);
2095  model->rho = Malloc(double,1);
2096  model->rho[0] = f.rho;
2097 
2098  int nSV = 0;
2099  int i;
2100  for(i=0;i<prob->l;i++)
2101  if(fabs(f.alpha[i]) > 0) ++nSV;
2102  model->l = nSV;
2103  model->SV = Malloc(svm_node *,nSV);
2104  model->sv_coef[0] = Malloc(double,nSV);
2105  int j = 0;
2106  for(i=0;i<prob->l;i++)
2107  if(fabs(f.alpha[i]) > 0)
2108  {
2109  model->SV[j] = prob->x[i];
2110  model->sv_coef[0][j] = f.alpha[i];
2111  ++j;
2112  }
2113 
2114  free(f.alpha);
2115  }
2116  else
2117  {
2118  // classification
2119  int l = prob->l;
2120  int nr_class;
2121  int *label = NULL;
2122  int *start = NULL;
2123  int *count = NULL;
2124  int *perm = Malloc(int,l);
2125 
2126  // group training data of the same class
2127  svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
2128  if(nr_class == 1)
2129  info("WARNING: training data in only one class. See README for details.\n");
2130 
2131  svm_node **x = Malloc(svm_node *,l);
2132  int i;
2133  for(i=0;i<l;i++)
2134  x[i] = prob->x[perm[i]];
2135 
2136  // calculate weighted C
2137 
2138  double *weighted_C = Malloc(double, nr_class);
2139  for(i=0;i<nr_class;i++)
2140  weighted_C[i] = param->C;
2141  for(i=0;i<param->nr_weight;i++)
2142  {
2143  int j;
2144  for(j=0;j<nr_class;j++)
2145  if(param->weight_label[i] == label[j])
2146  break;
2147  if(j == nr_class)
2148  fprintf(stderr,"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]);
2149  else
2150  weighted_C[j] *= param->weight[i];
2151  }
2152 
2153  // train k*(k-1)/2 models
2154 
2155  bool *nonzero = Malloc(bool,l);
2156  for(i=0;i<l;i++)
2157  nonzero[i] = false;
2158  decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
2159 
2160  double *probA=NULL,*probB=NULL;
2161  if (param->probability)
2162  {
2163  probA=Malloc(double,nr_class*(nr_class-1)/2);
2164  probB=Malloc(double,nr_class*(nr_class-1)/2);
2165  }
2166 
2167  int p = 0;
2168  for(i=0;i<nr_class;i++)
2169  for(int j=i+1;j<nr_class;j++)
2170  {
2171  svm_problem sub_prob;
2172  int si = start[i], sj = start[j];
2173  int ci = count[i], cj = count[j];
2174  sub_prob.l = ci+cj;
2175  sub_prob.x = Malloc(svm_node *,sub_prob.l);
2176  sub_prob.y = Malloc(double,sub_prob.l);
2177  int k;
2178  for(k=0;k<ci;k++)
2179  {
2180  sub_prob.x[k] = x[si+k];
2181  sub_prob.y[k] = +1;
2182  }
2183  for(k=0;k<cj;k++)
2184  {
2185  sub_prob.x[ci+k] = x[sj+k];
2186  sub_prob.y[ci+k] = -1;
2187  }
2188 
2189  if(param->probability)
2190  svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);
2191 
2192  f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
2193  for(k=0;k<ci;k++)
2194  if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
2195  nonzero[si+k] = true;
2196  for(k=0;k<cj;k++)
2197  if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
2198  nonzero[sj+k] = true;
2199  free(sub_prob.x);
2200  free(sub_prob.y);
2201  ++p;
2202  }
2203 
2204  // build output
2205 
2206  model->nr_class = nr_class;
2207 
2208  model->label = Malloc(int,nr_class);
2209  for(i=0;i<nr_class;i++)
2210  model->label[i] = label[i];
2211 
2212  model->rho = Malloc(double,nr_class*(nr_class-1)/2);
2213  for(i=0;i<nr_class*(nr_class-1)/2;i++)
2214  model->rho[i] = f[i].rho;
2215 
2216  if(param->probability)
2217  {
2218  model->probA = Malloc(double,nr_class*(nr_class-1)/2);
2219  model->probB = Malloc(double,nr_class*(nr_class-1)/2);
2220  for(i=0;i<nr_class*(nr_class-1)/2;i++)
2221  {
2222  model->probA[i] = probA[i];
2223  model->probB[i] = probB[i];
2224  }
2225  }
2226  else
2227  {
2228  model->probA=NULL;
2229  model->probB=NULL;
2230  }
2231 
2232  int total_sv = 0;
2233  int *nz_count = Malloc(int,nr_class);
2234  model->nSV = Malloc(int,nr_class);
2235  for(i=0;i<nr_class;i++)
2236  {
2237  int nSV = 0;
2238  for(int j=0;j<count[i];j++)
2239  if(nonzero[start[i]+j])
2240  {
2241  ++nSV;
2242  ++total_sv;
2243  }
2244  model->nSV[i] = nSV;
2245  nz_count[i] = nSV;
2246  }
2247 
2248  info("Total nSV = %d\n",total_sv);
2249 
2250  model->l = total_sv;
2251  model->SV = Malloc(svm_node *,total_sv);
2252  p = 0;
2253  for(i=0;i<l;i++)
2254  if(nonzero[i]) model->SV[p++] = x[i];
2255 
2256  int *nz_start = Malloc(int,nr_class);
2257  nz_start[0] = 0;
2258  for(i=1;i<nr_class;i++)
2259  nz_start[i] = nz_start[i-1]+nz_count[i-1];
2260 
2261  model->sv_coef = Malloc(double *,nr_class-1);
2262  for(i=0;i<nr_class-1;i++)
2263  model->sv_coef[i] = Malloc(double,total_sv);
2264 
2265  p = 0;
2266  for(i=0;i<nr_class;i++)
2267  for(int j=i+1;j<nr_class;j++)
2268  {
2269  // classifier (i,j): coefficients with
2270  // i are in sv_coef[j-1][nz_start[i]...],
2271  // j are in sv_coef[i][nz_start[j]...]
2272 
2273  int si = start[i];
2274  int sj = start[j];
2275  int ci = count[i];
2276  int cj = count[j];
2277 
2278  int q = nz_start[i];
2279  int k;
2280  for(k=0;k<ci;k++)
2281  if(nonzero[si+k])
2282  model->sv_coef[j-1][q++] = f[p].alpha[k];
2283  q = nz_start[j];
2284  for(k=0;k<cj;k++)
2285  if(nonzero[sj+k])
2286  model->sv_coef[i][q++] = f[p].alpha[ci+k];
2287  ++p;
2288  }
2289 
2290  free(label);
2291  free(probA);
2292  free(probB);
2293  free(count);
2294  free(perm);
2295  free(start);
2296  free(x);
2297  free(weighted_C);
2298  free(nonzero);
2299  for(i=0;i<nr_class*(nr_class-1)/2;i++)
2300  free(f[i].alpha);
2301  free(f);
2302  free(nz_count);
2303  free(nz_start);
2304  }
2305  return model;
2306 }
2307 
2308 // Stratified cross validation
2309 void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target)
2310 {
2311  int i;
2312  int *fold_start = Malloc(int,nr_fold+1);
2313  int l = prob->l;
2314  int *perm = Malloc(int,l);
2315  int nr_class;
2316 
2317  // stratified cv may not give leave-one-out rate
2318  // Each class to l folds -> some folds may have zero elements
2319  if((param->svm_type == C_SVC ||
2320  param->svm_type == NU_SVC) && nr_fold < l)
2321  {
2322  int *start = NULL;
2323  int *label = NULL;
2324  int *count = NULL;
2325  svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
2326 
2327  // random shuffle and then data grouped by fold using the array perm
2328  int *fold_count = Malloc(int,nr_fold);
2329  int c;
2330  int *index = Malloc(int,l);
2331  for(i=0;i<l;i++)
2332  index[i]=perm[i];
2333  for (c=0; c<nr_class; c++)
2334  for(i=0;i<count[c];i++)
2335  {
2336  int j = i+rand()%(count[c]-i);
2337  std::swap(index[start[c]+j],index[start[c]+i]);
2338  }
2339  for(i=0;i<nr_fold;i++)
2340  {
2341  fold_count[i] = 0;
2342  for (c=0; c<nr_class;c++)
2343  fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold;
2344  }
2345  fold_start[0]=0;
2346  for (i=1;i<=nr_fold;i++)
2347  fold_start[i] = fold_start[i-1]+fold_count[i-1];
2348  for (c=0; c<nr_class;c++)
2349  for(i=0;i<nr_fold;i++)
2350  {
2351  int begin = start[c]+i*count[c]/nr_fold;
2352  int end = start[c]+(i+1)*count[c]/nr_fold;
2353  for(int j=begin;j<end;j++)
2354  {
2355  perm[fold_start[i]] = index[j];
2356  fold_start[i]++;
2357  }
2358  }
2359  fold_start[0]=0;
2360  for (i=1;i<=nr_fold;i++)
2361  fold_start[i] = fold_start[i-1]+fold_count[i-1];
2362  free(start);
2363  free(label);
2364  free(count);
2365  free(index);
2366  free(fold_count);
2367  }
2368  else
2369  {
2370  for(i=0;i<l;i++) perm[i]=i;
2371  for(i=0;i<l;i++)
2372  {
2373  int j = i+rand()%(l-i);
2374  std::swap(perm[i],perm[j]);
2375  }
2376  for(i=0;i<=nr_fold;i++)
2377  fold_start[i]=i*l/nr_fold;
2378  }
2379 
2380  for(i=0;i<nr_fold;i++)
2381  {
2382  int begin = fold_start[i];
2383  int end = fold_start[i+1];
2384  int j,k;
2385  struct svm_problem subprob;
2386 
2387  subprob.l = l-(end-begin);
2388  subprob.x = Malloc(struct svm_node*,subprob.l);
2389  subprob.y = Malloc(double,subprob.l);
2390 
2391  k=0;
2392  for(j=0;j<begin;j++)
2393  {
2394  subprob.x[k] = prob->x[perm[j]];
2395  subprob.y[k] = prob->y[perm[j]];
2396  ++k;
2397  }
2398  for(j=end;j<l;j++)
2399  {
2400  subprob.x[k] = prob->x[perm[j]];
2401  subprob.y[k] = prob->y[perm[j]];
2402  ++k;
2403  }
2404  struct svm_model *submodel = svm_train(&subprob,param);
2405  if(param->probability &&
2406  (param->svm_type == C_SVC || param->svm_type == NU_SVC))
2407  {
2408  double *prob_estimates=Malloc(double,svm_get_nr_class(submodel));
2409  for(j=begin;j<end;j++)
2410  target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
2411  free(prob_estimates);
2412  }
2413  else
2414  for(j=begin;j<end;j++)
2415  target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
2416  svm_free_and_destroy_model(&submodel);
2417  free(subprob.x);
2418  free(subprob.y);
2419  }
2420  free(fold_start);
2421  free(perm);
2422 }
2423 
2424 
2425 int svm_get_svm_type(const svm_model *model)
2426 {
2427  return model->param.svm_type;
2428 }
2429 
2430 int svm_get_nr_class(const svm_model *model)
2431 {
2432  return model->nr_class;
2433 }
2434 
2435 void svm_get_labels(const svm_model *model, int* label)
2436 {
2437  if (model->label != NULL)
2438  for(int i=0;i<model->nr_class;i++)
2439  label[i] = model->label[i];
2440 }
2441 
2442 double svm_get_svr_probability(const svm_model *model)
2443 {
2444  if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
2445  model->probA!=NULL)
2446  return model->probA[0];
2447  else
2448  {
2449  fprintf(stderr,"Model doesn't contain information for SVR probability inference\n");
2450  return 0;
2451  }
2452 }
2453 
2454 double svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
2455 {
2456  int i;
2457  if(model->param.svm_type == ONE_CLASS ||
2458  model->param.svm_type == EPSILON_SVR ||
2459  model->param.svm_type == NU_SVR)
2460  {
2461  double *sv_coef = model->sv_coef[0];
2462  double sum = 0;
2463  for(i=0;i<model->l;i++)
2464  sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
2465  sum -= model->rho[0];
2466  *dec_values = sum;
2467 
2468  if(model->param.svm_type == ONE_CLASS)
2469  return (sum>0)?1:-1;
2470  else
2471  return sum;
2472  }
2473  else
2474  {
2475  int nr_class = model->nr_class;
2476  int l = model->l;
2477 
2478  double *kvalue = Malloc(double,l);
2479  for(i=0;i<l;i++)
2480  kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);
2481 
2482  int *start = Malloc(int,nr_class);
2483  start[0] = 0;
2484  for(i=1;i<nr_class;i++)
2485  start[i] = start[i-1]+model->nSV[i-1];
2486 
2487  int *vote = Malloc(int,nr_class);
2488  for(i=0;i<nr_class;i++)
2489  vote[i] = 0;
2490 
2491  int p=0;
2492  for(i=0;i<nr_class;i++)
2493  for(int j=i+1;j<nr_class;j++)
2494  {
2495  double sum = 0;
2496  int si = start[i];
2497  int sj = start[j];
2498  int ci = model->nSV[i];
2499  int cj = model->nSV[j];
2500 
2501  int k;
2502  double *coef1 = model->sv_coef[j-1];
2503  double *coef2 = model->sv_coef[i];
2504  for(k=0;k<ci;k++)
2505  sum += coef1[si+k] * kvalue[si+k];
2506  for(k=0;k<cj;k++)
2507  sum += coef2[sj+k] * kvalue[sj+k];
2508  sum -= model->rho[p];
2509  dec_values[p] = sum;
2510 
2511  if(dec_values[p] > 0)
2512  ++vote[i];
2513  else
2514  ++vote[j];
2515  p++;
2516  }
2517 
2518  int vote_max_idx = 0;
2519  for(i=1;i<nr_class;i++)
2520  if(vote[i] > vote[vote_max_idx])
2521  vote_max_idx = i;
2522 
2523  free(kvalue);
2524  free(start);
2525  free(vote);
2526  return model->label[vote_max_idx];
2527  }
2528 }
2529 
2530 double svm_predict(const svm_model *model, const svm_node *x)
2531 {
2532  int nr_class = model->nr_class;
2533  double *dec_values;
2534  if(model->param.svm_type == ONE_CLASS ||
2535  model->param.svm_type == EPSILON_SVR ||
2536  model->param.svm_type == NU_SVR)
2537  dec_values = Malloc(double, 1);
2538  else
2539  dec_values = Malloc(double, nr_class*(nr_class-1)/2);
2540  double pred_result = svm_predict_values(model, x, dec_values);
2541  free(dec_values);
2542  return pred_result;
2543 }
2544 
2545 double svm_predict_probability(
2546  const svm_model *model, const svm_node *x, double *prob_estimates)
2547 {
2548  if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
2549  model->probA!=NULL && model->probB!=NULL)
2550  {
2551  int i;
2552  int nr_class = model->nr_class;
2553  double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
2554  svm_predict_values(model, x, dec_values);
2555 
2556  double min_prob=1e-7;
2557  double **pairwise_prob=Malloc(double *,nr_class);
2558  for(i=0;i<nr_class;i++)
2559  pairwise_prob[i]=Malloc(double,nr_class);
2560  int k=0;
2561  for(i=0;i<nr_class;i++)
2562  for(int j=i+1;j<nr_class;j++)
2563  {
2564  pairwise_prob[i][j]=std::min(std::max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
2565  pairwise_prob[j][i]=1-pairwise_prob[i][j];
2566  k++;
2567  }
2568  multiclass_probability(nr_class,pairwise_prob,prob_estimates);
2569 
2570  int prob_max_idx = 0;
2571  for(i=1;i<nr_class;i++)
2572  if(prob_estimates[i] > prob_estimates[prob_max_idx])
2573  prob_max_idx = i;
2574  for(i=0;i<nr_class;i++)
2575  free(pairwise_prob[i]);
2576  free(dec_values);
2577  free(pairwise_prob);
2578  return model->label[prob_max_idx];
2579  }
2580  else
2581  return svm_predict(model, x);
2582 }
2583 
2584 static const char *svm_type_table[] =
2585 {
2586  "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
2587 };
2588 
2589 static const char *kernel_type_table[]=
2590 {
2591  "linear","polynomial","rbf","sigmoid","precomputed",NULL
2592 };
2593 
2594 int svm_save_model(const char *model_file_name, const svm_model *model)
2595 {
2596  FILE *fp = fopen(model_file_name,"w");
2597  if(fp==NULL) return -1;
2598 
2599  char *old_locale = strdup(setlocale(LC_ALL, NULL));
2600  setlocale(LC_ALL, "C");
2601 
2602  const svm_parameter& param = model->param;
2603 
2604  fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
2605  fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);
2606 
2607  if(param.kernel_type == POLY)
2608  fprintf(fp,"degree %d\n", param.degree);
2609 
2610  if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
2611  fprintf(fp,"gamma %g\n", param.gamma);
2612 
2613  if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
2614  fprintf(fp,"coef0 %g\n", param.coef0);
2615 
2616  int nr_class = model->nr_class;
2617  int l = model->l;
2618  fprintf(fp, "nr_class %d\n", nr_class);
2619  fprintf(fp, "total_sv %d\n",l);
2620 
2621  {
2622  fprintf(fp, "rho");
2623  for(int i=0;i<nr_class*(nr_class-1)/2;i++)
2624  fprintf(fp," %g",model->rho[i]);
2625  fprintf(fp, "\n");
2626  }
2627 
2628  if(model->label)
2629  {
2630  fprintf(fp, "label");
2631  for(int i=0;i<nr_class;i++)
2632  fprintf(fp," %d",model->label[i]);
2633  fprintf(fp, "\n");
2634  }
2635 
2636  if(model->probA) // regression has probA only
2637  {
2638  fprintf(fp, "probA");
2639  for(int i=0;i<nr_class*(nr_class-1)/2;i++)
2640  fprintf(fp," %g",model->probA[i]);
2641  fprintf(fp, "\n");
2642  }
2643  if(model->probB)
2644  {
2645  fprintf(fp, "probB");
2646  for(int i=0;i<nr_class*(nr_class-1)/2;i++)
2647  fprintf(fp," %g",model->probB[i]);
2648  fprintf(fp, "\n");
2649  }
2650 
2651  if(model->nSV)
2652  {
2653  fprintf(fp, "nr_sv");
2654  for(int i=0;i<nr_class;i++)
2655  fprintf(fp," %d",model->nSV[i]);
2656  fprintf(fp, "\n");
2657  }
2658 
2659  fprintf(fp, "SV\n");
2660  const double * const *sv_coef = model->sv_coef;
2661  const svm_node * const *SV = model->SV;
2662 
2663  for(int i=0;i<l;i++)
2664  {
2665  for(int j=0;j<nr_class-1;j++)
2666  fprintf(fp, "%.16g ",sv_coef[j][i]);
2667 
2668  const svm_node *p = SV[i];
2669 
2670  if(param.kernel_type == PRECOMPUTED)
2671  fprintf(fp,"0:%d ",(int)(p->value));
2672  else
2673  while(p->index != -1)
2674  {
2675  fprintf(fp,"%d:%.8g ",p->index,p->value);
2676  p++;
2677  }
2678  fprintf(fp, "\n");
2679  }
2680 
2681  setlocale(LC_ALL, old_locale);
2682  free(old_locale);
2683 
2684  if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
2685  else return 0;
2686 }
2687 
2688 static char *line = NULL;
2689 static int max_line_len;
2690 
2691 static char* readline(FILE *input)
2692 {
2693  int len;
2694 
2695  if(fgets(line,max_line_len,input) == NULL)
2696  return NULL;
2697 
2698  while(strrchr(line,'\n') == NULL)
2699  {
2700  max_line_len *= 2;
2701  line = (char *) realloc(line,max_line_len);
2702  len = (int) strlen(line);
2703  if(fgets(line+len,max_line_len-len,input) == NULL)
2704  break;
2705  }
2706  return line;
2707 }
2708 
2709 svm_model *svm_load_model(const char *model_file_name)
2710 {
2711  FILE *fp = fopen(model_file_name,"rb");
2712  if(fp==NULL) return NULL;
2713 
2714  char *old_locale = strdup(setlocale(LC_ALL, NULL));
2715  setlocale(LC_ALL, "C");
2716 
2717  // read parameters
2718 
2719  svm_model *model = Malloc(svm_model,1);
2720  svm_parameter& param = model->param;
2721  model->rho = NULL;
2722  model->probA = NULL;
2723  model->probB = NULL;
2724  model->label = NULL;
2725  model->nSV = NULL;
2726 
2727  char cmd[81];
2728  for(unsigned int i=0; i<81; i++) cmd[i] = '\0';
2729  while(1)
2730  {
2731  fscanf(fp,"%80s",cmd);
2732 
2733  if(strcmp(cmd,"svm_type")==0)
2734  {
2735  fscanf(fp,"%80s",cmd);
2736  int i;
2737  for(i=0;svm_type_table[i];i++)
2738  {
2739  if(strcmp(svm_type_table[i],cmd)==0)
2740  {
2741  param.svm_type=i;
2742  break;
2743  }
2744  }
2745  if(svm_type_table[i] == NULL)
2746  {
2747  fprintf(stderr,"unknown svm type.\n");
2748 
2749  setlocale(LC_ALL, old_locale);
2750  free(old_locale);
2751  free(model->rho);
2752  free(model->label);
2753  free(model->nSV);
2754  free(model);
2755  return NULL;
2756  }
2757  }
2758  else if(strcmp(cmd,"kernel_type")==0)
2759  {
2760  fscanf(fp,"%80s",cmd);
2761  int i;
2762  for(i=0;kernel_type_table[i];i++)
2763  {
2764  if(strcmp(kernel_type_table[i],cmd)==0)
2765  {
2766  param.kernel_type=i;
2767  break;
2768  }
2769  }
2770  if(kernel_type_table[i] == NULL)
2771  {
2772  fprintf(stderr,"unknown kernel function.\n");
2773 
2774  setlocale(LC_ALL, old_locale);
2775  free(old_locale);
2776  free(model->rho);
2777  free(model->label);
2778  free(model->nSV);
2779  free(model);
2780  return NULL;
2781  }
2782  }
2783  else if(strcmp(cmd,"degree")==0)
2784  fscanf(fp,"%d",&param.degree);
2785  else if(strcmp(cmd,"gamma")==0)
2786  fscanf(fp,"%lf",&param.gamma);
2787  else if(strcmp(cmd,"coef0")==0)
2788  fscanf(fp,"%lf",&param.coef0);
2789  else if(strcmp(cmd,"nr_class")==0)
2790  fscanf(fp,"%d",&model->nr_class);
2791  else if(strcmp(cmd,"total_sv")==0)
2792  fscanf(fp,"%d",&model->l);
2793  else if(strcmp(cmd,"rho")==0)
2794  {
2795  int n = model->nr_class * (model->nr_class-1)/2;
2796  model->rho = Malloc(double,n);
2797  for(int i=0;i<n;i++)
2798  fscanf(fp,"%lf",&model->rho[i]);
2799  }
2800  else if(strcmp(cmd,"label")==0)
2801  {
2802  int n = model->nr_class;
2803  model->label = Malloc(int,n);
2804  for(int i=0;i<n;i++)
2805  fscanf(fp,"%d",&model->label[i]);
2806  }
2807  else if(strcmp(cmd,"probA")==0)
2808  {
2809  int n = model->nr_class * (model->nr_class-1)/2;
2810  model->probA = Malloc(double,n);
2811  for(int i=0;i<n;i++)
2812  fscanf(fp,"%lf",&model->probA[i]);
2813  }
2814  else if(strcmp(cmd,"probB")==0)
2815  {
2816  int n = model->nr_class * (model->nr_class-1)/2;
2817  model->probB = Malloc(double,n);
2818  for(int i=0;i<n;i++)
2819  fscanf(fp,"%lf",&model->probB[i]);
2820  }
2821  else if(strcmp(cmd,"nr_sv")==0)
2822  {
2823  int n = model->nr_class;
2824  model->nSV = Malloc(int,n);
2825  for(int i=0;i<n;i++)
2826  fscanf(fp,"%d",&model->nSV[i]);
2827  }
2828  else if(strcmp(cmd,"SV")==0)
2829  {
2830  while(1)
2831  {
2832  int c = getc(fp);
2833  if(c==EOF || c=='\n') break;
2834  }
2835  break;
2836  }
2837  else
2838  {
2839  fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
2840 
2841  setlocale(LC_ALL, old_locale);
2842  free(old_locale);
2843  free(model->rho);
2844  free(model->label);
2845  free(model->nSV);
2846  free(model);
2847  return NULL;
2848  }
2849  }
2850 
2851  // read sv_coef and SV
2852 
2853  int elements = 0;
2854  long pos = ftell(fp);
2855 
2856  max_line_len = 1024;
2857  line = Malloc(char,max_line_len);
2858  char *p,*endptr,*idx,*val;
2859 
2860  while(readline(fp)!=NULL)
2861  {
2862  p = strtok(line,":");
2863  while(1)
2864  {
2865  p = strtok(NULL,":");
2866  if(p == NULL)
2867  break;
2868  ++elements;
2869  }
2870  }
2871  elements += model->l;
2872 
2873  fseek(fp,pos,SEEK_SET);
2874 
2875  int m = model->nr_class - 1;
2876  int l = model->l;
2877  model->sv_coef = Malloc(double *,m);
2878  int i;
2879  for(i=0;i<m;i++)
2880  model->sv_coef[i] = Malloc(double,l);
2881  model->SV = Malloc(svm_node*,l);
2882  svm_node *x_space = NULL;
2883  if(l>0) x_space = Malloc(svm_node,elements);
2884 
2885  int j=0;
2886  for(i=0;i<l;i++)
2887  {
2888  readline(fp);
2889  model->SV[i] = &x_space[j];
2890 
2891  p = strtok(line, " \t");
2892  model->sv_coef[0][i] = strtod(p,&endptr);
2893  for(int k=1;k<m;k++)
2894  {
2895  p = strtok(NULL, " \t");
2896  model->sv_coef[k][i] = strtod(p,&endptr);
2897  }
2898 
2899  while(1)
2900  {
2901  idx = strtok(NULL, ":");
2902  val = strtok(NULL, " \t");
2903 
2904  if(val == NULL)
2905  break;
2906  x_space[j].index = (int) strtol(idx,&endptr,10);
2907  x_space[j].value = strtod(val,&endptr);
2908 
2909  ++j;
2910  }
2911  x_space[j++].index = -1;
2912  }
2913  free(line);
2914 
2915  setlocale(LC_ALL, old_locale);
2916  free(old_locale);
2917 
2918  if (ferror(fp) != 0 || fclose(fp) != 0)
2919  return NULL;
2920 
2921  model->free_sv = 1; // XXX
2922  return model;
2923 }
2924 
2925 void svm_free_model_content(svm_model* model_ptr)
2926 {
2927  if(model_ptr->free_sv && model_ptr->l > 0 && model_ptr->SV != NULL)
2928  free((void *)(model_ptr->SV[0]));
2929  if(model_ptr->sv_coef)
2930  {
2931  for(int i=0;i<model_ptr->nr_class-1;i++)
2932  free(model_ptr->sv_coef[i]);
2933  }
2934 
2935  free(model_ptr->SV);
2936  model_ptr->SV = NULL;
2937 
2938  free(model_ptr->sv_coef);
2939  model_ptr->sv_coef = NULL;
2940 
2941  free(model_ptr->rho);
2942  model_ptr->rho = NULL;
2943 
2944  free(model_ptr->label);
2945  model_ptr->label= NULL;
2946 
2947  free(model_ptr->probA);
2948  model_ptr->probA = NULL;
2949 
2950  free(model_ptr->probB);
2951  model_ptr->probB= NULL;
2952 
2953  free(model_ptr->nSV);
2954  model_ptr->nSV = NULL;
2955 }
2956 
2957 void svm_free_and_destroy_model(svm_model** model_ptr_ptr)
2958 {
2959  if(model_ptr_ptr != NULL && *model_ptr_ptr != NULL)
2960  {
2961  svm_free_model_content(*model_ptr_ptr);
2962  free(*model_ptr_ptr);
2963  *model_ptr_ptr = NULL;
2964  }
2965 }
2966 
2967 void svm_destroy_param(svm_parameter* param)
2968 {
2969  free(param->weight_label);
2970  free(param->weight);
2971 }
2972 
2973 const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
2974 {
2975  // svm_type
2976  int svm_type = param->svm_type;
2977  if(svm_type != C_SVC &&
2978  svm_type != NU_SVC &&
2979  svm_type != ONE_CLASS &&
2980  svm_type != EPSILON_SVR &&
2981  svm_type != NU_SVR)
2982  return "unknown svm type";
2983 
2984  // kernel_type, degree
2985  int kernel_type = param->kernel_type;
2986  if(kernel_type != LINEAR &&
2987  kernel_type != POLY &&
2988  kernel_type != RBF &&
2989  kernel_type != SIGMOID &&
2990  kernel_type != PRECOMPUTED)
2991  return "unknown kernel type";
2992 
2993  if(param->gamma < 0)
2994  return "gamma < 0";
2995 
2996  if(param->degree < 0)
2997  return "degree of polynomial kernel < 0";
2998 
2999  // cache_size,eps,C,nu,p,shrinking
3000 
3001  if(param->cache_size <= 0)
3002  return "cache_size <= 0";
3003 
3004  if(param->eps <= 0)
3005  return "eps <= 0";
3006 
3007  if(svm_type == C_SVC ||
3008  svm_type == EPSILON_SVR ||
3009  svm_type == NU_SVR)
3010  if(param->C <= 0)
3011  return "C <= 0";
3012 
3013  if(svm_type == NU_SVC ||
3014  svm_type == ONE_CLASS ||
3015  svm_type == NU_SVR)
3016  if(param->nu <= 0 || param->nu > 1)
3017  return "nu <= 0 or nu > 1";
3018 
3019  if(svm_type == EPSILON_SVR)
3020  if(param->p < 0)
3021  return "p < 0";
3022 
3023  if(param->shrinking != 0 &&
3024  param->shrinking != 1)
3025  return "shrinking != 0 and shrinking != 1";
3026 
3027  if(param->probability != 0 &&
3028  param->probability != 1)
3029  return "probability != 0 and probability != 1";
3030 
3031  if(param->probability == 1 &&
3032  svm_type == ONE_CLASS)
3033  return "one-class SVM probability output not supported yet";
3034 
3035 
3036  // check whether nu-svc is feasible
3037  if(svm_type == NU_SVC)
3038  {
3039  int l = prob->l;
3040  int max_nr_class = 16;
3041  int nr_class = 0;
3042  int *label = Malloc(int,max_nr_class);
3043  int *count = Malloc(int,max_nr_class);
3044 
3045  int i;
3046  for(i=0;i<l;i++)
3047  {
3048  int this_label = (int)prob->y[i];
3049  int j;
3050  for(j=0;j<nr_class;j++)
3051  if(this_label == label[j])
3052  {
3053  ++count[j];
3054  break;
3055  }
3056  if(j == nr_class)
3057  {
3058  if(nr_class == max_nr_class)
3059  {
3060  max_nr_class *= 2;
3061  label = (int *)realloc(label,max_nr_class*sizeof(int));
3062  count = (int *)realloc(count,max_nr_class*sizeof(int));
3063  }
3064  label[nr_class] = this_label;
3065  count[nr_class] = 1;
3066  ++nr_class;
3067  }
3068  }
3069 
3070  for(i=0;i<nr_class;i++)
3071  {
3072  int n1 = count[i];
3073  for(int j=i+1;j<nr_class;j++)
3074  {
3075  int n2 = count[j];
3076  if(param->nu*(n1+n2)/2 > std::min(n1,n2))
3077  {
3078  free(label);
3079  free(count);
3080  return "specified nu is infeasible";
3081  }
3082  }
3083  }
3084  free(label);
3085  free(count);
3086  }
3087 
3088  return NULL;
3089 }
3090 
3091 int svm_check_probability_model(const svm_model *model)
3092 {
3093  return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
3094  model->probA!=NULL && model->probB!=NULL) ||
3095  ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
3096  model->probA!=NULL);
3097 }
3098 
3099 void svm_set_print_string_function(void (*print_func)(const char *))
3100 {
3101  if(print_func == NULL)
3102  svm_print_string = &print_string_stdout;
3103  else
3104  svm_print_string = print_func;
3105 }
3106 
3107 } //End of namespace LIBSVM
Definition: libsvm.cpp:4