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