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