SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
libocas.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * libocas.c: Implementation of the OCAS solver for training
8  * linear SVM classifiers.
9  *
10  * Copyright (C) 2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
11  * Soeren Sonnenburg, soeren.sonnenburg@first.fraunhofer.de
12  *-------------------------------------------------------------------- */
13 
14 #include <stdlib.h>
15 #include <string.h>
16 #include <math.h>
17 #include <sys/time.h>
18 #include <time.h>
19 #include <stdio.h>
20 #include <stdint.h>
21 
25 
26 namespace shogun
27 {
28 #define MU 0.1 /* must be from (0,1> 1..means that OCAS becomes equivalent to CPA */
29  /* see paper Franc&Sonneburg JMLR 2009 */
30 
31 static const uint32_t QPSolverMaxIter = 10000000;
32 
33 static float64_t *H;
34 static uint32_t BufSize;
35 
36 /*----------------------------------------------------------------------
37  Returns pointer at i-th column of Hessian matrix.
38  ----------------------------------------------------------------------*/
39 static const float64_t *get_col( uint32_t i)
40 {
41  return( &H[ BufSize*i ] );
42 }
43 
44 /*----------------------------------------------------------------------
45  Returns time of the day in seconds.
46  ----------------------------------------------------------------------*/
48 {
49  struct timeval tv;
50  if (gettimeofday(&tv, NULL)==0)
51  return tv.tv_sec+((float64_t)(tv.tv_usec))/1e6;
52  else
53  return 0.0;
54 }
55 
56 
57 /*----------------------------------------------------------------------
58  Linear binary Ocas-SVM solver with additinal contraint enforceing
59  a subset of weights (indices of the weights given by num_nnw/nnw_idx)
60  to be non-negative.
61 
62  ----------------------------------------------------------------------*/
63 ocas_return_value_T svm_ocas_solver_nnw(
64  float64_t C,
65  uint32_t nData,
66  uint32_t num_nnw,
67  uint32_t* nnw_idx,
68  float64_t TolRel,
69  float64_t TolAbs,
70  float64_t QPBound,
71  float64_t MaxTime,
72  uint32_t _BufSize,
73  uint8_t Method,
74  int (*add_pw_constr)(uint32_t, uint32_t, void*),
75  void (*clip_neg_W)(uint32_t, uint32_t*, void*),
76  void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
77  float64_t (*update_W)(float64_t, void*),
78  int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*),
79  int (*compute_output)(float64_t*, void* ),
80  int (*sort)(float64_t*, float64_t*, uint32_t),
81  void (*ocas_print)(ocas_return_value_T),
82  void* user_data)
83 {
84  ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
85  float64_t *b, *alpha, *diag_H;
86  float64_t *output, *old_output;
87  float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
88  float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
89  float64_t start_time, ocas_start_time;
90  uint32_t cut_length;
91  uint32_t i, *new_cut;
92  uint32_t *I;
93 /* uint8_t S = 1; */
94  libqp_state_T qp_exitflag;
95 
96  float64_t max_cp_norm;
97  float64_t max_b;
98  float64_t _C[2];
99  uint8_t S[2];
100 
101  ocas_start_time = get_time();
102  ocas.qp_solver_time = 0;
103  ocas.output_time = 0;
104  ocas.sort_time = 0;
105  ocas.add_time = 0;
106  ocas.w_time = 0;
107  ocas.print_time = 0;
108 
109  BufSize = _BufSize;
110 
111  QPSolverTolRel = TolRel*0.5;
112 
113  H=NULL;
114  b=NULL;
115  alpha=NULL;
116  new_cut=NULL;
117  I=NULL;
118  diag_H=NULL;
119  output=NULL;
120  old_output=NULL;
121  hpf=NULL;
122  hpb = NULL;
123  Ci=NULL;
124  Bi=NULL;
125 
126  /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
127  H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
128  if(H == NULL)
129  {
130  ocas.exitflag=-2;
131  goto cleanup;
132  }
133 
134  /* bias of cutting planes */
135  b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
136  if(b == NULL)
137  {
138  ocas.exitflag=-2;
139  goto cleanup;
140  }
141 
142  alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
143  if(alpha == NULL)
144  {
145  ocas.exitflag=-2;
146  goto cleanup;
147  }
148 
149  /* indices of examples which define a new cut */
150  new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
151  if(new_cut == NULL)
152  {
153  ocas.exitflag=-2;
154  goto cleanup;
155  }
156 
157  I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
158  if(I == NULL)
159  {
160  ocas.exitflag=-2;
161  goto cleanup;
162  }
163 
164  for(i=0; i< BufSize; i++) I[i] = 2;
165 
166  diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
167  if(diag_H == NULL)
168  {
169  ocas.exitflag=-2;
170  goto cleanup;
171  }
172 
173  output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
174  if(output == NULL)
175  {
176  ocas.exitflag=-2;
177  goto cleanup;
178  }
179 
180  old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
181  if(old_output == NULL)
182  {
183  ocas.exitflag=-2;
184  goto cleanup;
185  }
186 
187  /* array of hinge points used in line-serach */
188  hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
189  if(hpf == NULL)
190  {
191  ocas.exitflag=-2;
192  goto cleanup;
193  }
194 
195  hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
196  if(hpb == NULL)
197  {
198  ocas.exitflag=-2;
199  goto cleanup;
200  }
201 
202  /* vectors Ci, Bi are used in the line search procedure */
203  Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
204  if(Ci == NULL)
205  {
206  ocas.exitflag=-2;
207  goto cleanup;
208  }
209 
210  Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
211  if(Bi == NULL)
212  {
213  ocas.exitflag=-2;
214  goto cleanup;
215  }
216 
217  /* initial cutting planes implementing the non-negativity constraints on W*/
218  _C[0]=10000000.0;
219  _C[1]=C;
220  S[0]=1;
221  S[1]=1;
222  for(i=0; i < num_nnw; i++)
223  {
224  if(add_pw_constr(nnw_idx[i],i, user_data) != 0)
225  {
226  ocas.exitflag=-2;
227  goto cleanup;
228  }
229  diag_H[i] = 1.0;
230  H[LIBOCAS_INDEX(i,i,BufSize)] = 1.0;
231  I[i] = 1;
232  }
233 
234  max_cp_norm = 1;
235  max_b = 0;
236 
237  /* */
238  ocas.nCutPlanes = num_nnw;
239  ocas.exitflag = 0;
240  ocas.nIter = 0;
241 
242  /* Compute initial value of Q_P assuming that W is zero vector.*/
243  sq_norm_W = 0;
244  xi = nData;
245  ocas.Q_P = 0.5*sq_norm_W + C*xi;
246  ocas.Q_D = 0;
247 
248  /* Compute the initial cutting plane */
249  cut_length = nData;
250  for(i=0; i < nData; i++)
251  new_cut[i] = i;
252 
253  ocas.trn_err = nData;
254  ocas.ocas_time = get_time() - ocas_start_time;
255  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n",
256  ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P));
257  */
258  ocas_print(ocas);
259 
260  /* main loop */
261  while( ocas.exitflag == 0 )
262  {
263  ocas.nIter++;
264 
265  /* append a new cut to the buffer and update H */
266  b[ocas.nCutPlanes] = -(float64_t)cut_length;
267 
268  max_b = LIBOCAS_MAX(max_b,(float64_t)cut_length);
269 
270  start_time = get_time();
271 
272  if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
273  {
274  ocas.exitflag=-2;
275  goto cleanup;
276  }
277 
278  ocas.add_time += get_time() - start_time;
279 
280  /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
281  diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
282  for(i=0; i < ocas.nCutPlanes; i++) {
283  H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
284  }
285 
286  max_cp_norm = LIBOCAS_MAX(max_cp_norm, sqrt(diag_H[ocas.nCutPlanes]));
287 
288 
289  ocas.nCutPlanes++;
290 
291  /* call inner QP solver */
292  start_time = get_time();
293 
294  /* compute upper bound on sum of dual variables associated with the positivity constraints */
295  _C[0] = sqrt((float64_t)nData)*(sqrt(C)*sqrt(max_b) + C*max_cp_norm);
296 
297 /* qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
298  ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);*/
299  qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, _C, I, S, alpha,
300  ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
301 
302  ocas.qp_exitflag = qp_exitflag.exitflag;
303 
304  ocas.qp_solver_time += get_time() - start_time;
305  ocas.Q_D = -qp_exitflag.QP;
306 
307  ocas.nNZAlpha = 0;
308  for(i=0; i < ocas.nCutPlanes; i++) {
309  if( alpha[i] != 0) ocas.nNZAlpha++;
310  }
311 
312  sq_norm_oldW = sq_norm_W;
313  start_time = get_time();
314  compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
315  clip_neg_W(num_nnw, nnw_idx, user_data);
316  ocas.w_time += get_time() - start_time;
317 
318  /* select a new cut */
319  switch( Method )
320  {
321  /* cutting plane algorithm implemented in SVMperf and BMRM */
322  case 0:
323 
324  start_time = get_time();
325  if( compute_output( output, user_data ) != 0)
326  {
327  ocas.exitflag=-2;
328  goto cleanup;
329  }
330  ocas.output_time += get_time()-start_time;
331 
332  xi = 0;
333  cut_length = 0;
334  ocas.trn_err = 0;
335  for(i=0; i < nData; i++)
336  {
337  if(output[i] <= 0) ocas.trn_err++;
338 
339  if(output[i] <= 1) {
340  xi += 1 - output[i];
341  new_cut[cut_length] = i;
342  cut_length++;
343  }
344  }
345  ocas.Q_P = 0.5*sq_norm_W + C*xi;
346 
347  ocas.ocas_time = get_time() - ocas_start_time;
348 
349  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
350  ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
351  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
352  */
353 
354  start_time = get_time();
355  ocas_print(ocas);
356  ocas.print_time += get_time() - start_time;
357 
358  break;
359 
360 
361  /* Ocas strategy */
362  case 1:
363 
364  /* Linesearch */
365  A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
366  B0 = dot_prod_WoldW - sq_norm_oldW;
367 
368  memcpy( old_output, output, sizeof(float64_t)*nData );
369 
370  start_time = get_time();
371  if( compute_output( output, user_data ) != 0)
372  {
373  ocas.exitflag=-2;
374  goto cleanup;
375  }
376  ocas.output_time += get_time()-start_time;
377 
378  uint32_t num_hp = 0;
379  GradVal = B0;
380  for(i=0; i< nData; i++) {
381 
382  Ci[i] = C*(1-old_output[i]);
383  Bi[i] = C*(old_output[i] - output[i]);
384 
385  float64_t val;
386  if(Bi[i] != 0)
387  val = -Ci[i]/Bi[i];
388  else
389  val = -LIBOCAS_PLUS_INF;
390 
391  if (val>0)
392  {
393 /* hpi[num_hp] = i;*/
394  hpb[num_hp] = Bi[i];
395  hpf[num_hp] = val;
396  num_hp++;
397  }
398 
399  if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
400  GradVal += Bi[i];
401 
402  }
403 
404  t = 0;
405  if( GradVal < 0 )
406  {
407  start_time = get_time();
408 /* if( sort(hpf, hpi, num_hp) != 0)*/
409  if( sort(hpf, hpb, num_hp) != 0 )
410  {
411  ocas.exitflag=-2;
412  goto cleanup;
413  }
414  ocas.sort_time += get_time() - start_time;
415 
416  float64_t t_new, GradVal_new;
417  i = 0;
418  while( GradVal < 0 && i < num_hp )
419  {
420  t_new = hpf[i];
421  GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
422 
423  if( GradVal_new >= 0 )
424  {
425  t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
426  }
427  else
428  {
429  t = t_new;
430  i++;
431  }
432 
433  GradVal = GradVal_new;
434  }
435  }
436 
437  /*
438  t = hpf[0] - 1;
439  i = 0;
440  GradVal = t*A0 + Bsum;
441  while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) {
442  t = hpf[i];
443  Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]);
444  GradVal = t*A0 + Bsum;
445  i++;
446  }
447  */
448  t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
449 
450  /* this guarantees that the new solution will not violate the positivity constraints on W */
451  t = LIBOCAS_MIN(t,1);
452 
453  t1 = t; /* new (best so far) W */
454  t2 = t+MU*(1.0-t); /* new cutting plane */
455  /* t2 = t+(1.0-t)/10.0; */
456 
457  /* update W to be the best so far solution */
458  sq_norm_W = update_W( t1, user_data );
459 
460  /* select a new cut */
461  xi = 0;
462  cut_length = 0;
463  ocas.trn_err = 0;
464  for(i=0; i < nData; i++ ) {
465 
466  if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 )
467  {
468  new_cut[cut_length] = i;
469  cut_length++;
470  }
471 
472  output[i] = old_output[i]*(1-t1) + t1*output[i];
473 
474  if( output[i] <= 1) xi += 1-output[i];
475  if( output[i] <= 0) ocas.trn_err++;
476 
477  }
478 
479  ocas.Q_P = 0.5*sq_norm_W + C*xi;
480 
481  ocas.ocas_time = get_time() - ocas_start_time;
482 
483  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
484  ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
485  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
486  */
487 
488  start_time = get_time();
489  ocas_print(ocas);
490  ocas.print_time += get_time() - start_time;
491 
492  break;
493  }
494 
495  /* Stopping conditions */
496  if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
497  if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
498  if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
499  if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
500  if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
501 
502  } /* end of the main loop */
503 
504 cleanup:
505 
506  LIBOCAS_FREE(H);
507  LIBOCAS_FREE(b);
508  LIBOCAS_FREE(alpha);
509  LIBOCAS_FREE(new_cut);
510  LIBOCAS_FREE(I);
511  LIBOCAS_FREE(diag_H);
512  LIBOCAS_FREE(output);
513  LIBOCAS_FREE(old_output);
514  LIBOCAS_FREE(hpf);
515 /* LIBOCAS_FREE(hpi);*/
516  LIBOCAS_FREE(hpb);
517  LIBOCAS_FREE(Ci);
518  LIBOCAS_FREE(Bi);
519 
520  ocas.ocas_time = get_time() - ocas_start_time;
521 
522  return(ocas);
523 }
524 
525 
526 
527 /*----------------------------------------------------------------------
528  Linear binary Ocas-SVM solver.
529  ----------------------------------------------------------------------*/
530 ocas_return_value_T svm_ocas_solver(
531  float64_t C,
532  uint32_t nData,
533  float64_t TolRel,
534  float64_t TolAbs,
535  float64_t QPBound,
536  float64_t MaxTime,
537  uint32_t _BufSize,
538  uint8_t Method,
539  void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
540  float64_t (*update_W)(float64_t, void*),
541  int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*),
542  int (*compute_output)(float64_t*, void* ),
543  int (*sort)(float64_t*, float64_t*, uint32_t),
544  void (*ocas_print)(ocas_return_value_T),
545  void* user_data)
546 {
547  ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
548  float64_t *b, *alpha, *diag_H;
549  float64_t *output, *old_output;
550  float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
551  float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
552  float64_t start_time, ocas_start_time;
553  uint32_t cut_length;
554  uint32_t i, *new_cut;
555  uint32_t *I;
556  uint8_t S = 1;
557  libqp_state_T qp_exitflag;
558 
559  ocas_start_time = get_time();
560  ocas.qp_solver_time = 0;
561  ocas.output_time = 0;
562  ocas.sort_time = 0;
563  ocas.add_time = 0;
564  ocas.w_time = 0;
565  ocas.print_time = 0;
566  float64_t gap;
567 
568  BufSize = _BufSize;
569 
570  QPSolverTolRel = TolRel*0.5;
571 
572  H=NULL;
573  b=NULL;
574  alpha=NULL;
575  new_cut=NULL;
576  I=NULL;
577  diag_H=NULL;
578  output=NULL;
579  old_output=NULL;
580  hpf=NULL;
581  hpb = NULL;
582  Ci=NULL;
583  Bi=NULL;
584 
585  /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
586  H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
587  if(H == NULL)
588  {
589  ocas.exitflag=-2;
590  goto cleanup;
591  }
592 
593  /* bias of cutting planes */
594  b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
595  if(b == NULL)
596  {
597  ocas.exitflag=-2;
598  goto cleanup;
599  }
600 
601  alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
602  if(alpha == NULL)
603  {
604  ocas.exitflag=-2;
605  goto cleanup;
606  }
607 
608  /* indices of examples which define a new cut */
609  new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
610  if(new_cut == NULL)
611  {
612  ocas.exitflag=-2;
613  goto cleanup;
614  }
615 
616  I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
617  if(I == NULL)
618  {
619  ocas.exitflag=-2;
620  goto cleanup;
621  }
622 
623  for(i=0; i< BufSize; i++) I[i] = 1;
624 
625  diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
626  if(diag_H == NULL)
627  {
628  ocas.exitflag=-2;
629  goto cleanup;
630  }
631 
632  output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
633  if(output == NULL)
634  {
635  ocas.exitflag=-2;
636  goto cleanup;
637  }
638 
639  old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
640  if(old_output == NULL)
641  {
642  ocas.exitflag=-2;
643  goto cleanup;
644  }
645 
646  /* array of hinge points used in line-serach */
647  hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
648  if(hpf == NULL)
649  {
650  ocas.exitflag=-2;
651  goto cleanup;
652  }
653 
654  hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
655  if(hpb == NULL)
656  {
657  ocas.exitflag=-2;
658  goto cleanup;
659  }
660 
661  /* vectors Ci, Bi are used in the line search procedure */
662  Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
663  if(Ci == NULL)
664  {
665  ocas.exitflag=-2;
666  goto cleanup;
667  }
668 
669  Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
670  if(Bi == NULL)
671  {
672  ocas.exitflag=-2;
673  goto cleanup;
674  }
675 
676  ocas.nCutPlanes = 0;
677  ocas.exitflag = 0;
678  ocas.nIter = 0;
679 
680  /* Compute initial value of Q_P assuming that W is zero vector.*/
681  sq_norm_W = 0;
682  xi = nData;
683  ocas.Q_P = 0.5*sq_norm_W + C*xi;
684  ocas.Q_D = 0;
685 
686  /* Compute the initial cutting plane */
687  cut_length = nData;
688  for(i=0; i < nData; i++)
689  new_cut[i] = i;
690 
691  gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P);
692  SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6);
693 
694  ocas.trn_err = nData;
695  ocas.ocas_time = get_time() - ocas_start_time;
696  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n",
697  ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P));
698  */
699  ocas_print(ocas);
700 
701  /* main loop */
702  while( ocas.exitflag == 0 )
703  {
704  ocas.nIter++;
705 
706  /* append a new cut to the buffer and update H */
707  b[ocas.nCutPlanes] = -(float64_t)cut_length;
708 
709  start_time = get_time();
710 
711  if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
712  {
713  ocas.exitflag=-2;
714  goto cleanup;
715  }
716 
717  ocas.add_time += get_time() - start_time;
718 
719  /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
720  diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
721  for(i=0; i < ocas.nCutPlanes; i++) {
722  H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
723  }
724 
725  ocas.nCutPlanes++;
726 
727  /* call inner QP solver */
728  start_time = get_time();
729 
730  qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
731  ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
732 
733  ocas.qp_exitflag = qp_exitflag.exitflag;
734 
735  ocas.qp_solver_time += get_time() - start_time;
736  ocas.Q_D = -qp_exitflag.QP;
737 
738  ocas.nNZAlpha = 0;
739  for(i=0; i < ocas.nCutPlanes; i++) {
740  if( alpha[i] != 0) ocas.nNZAlpha++;
741  }
742 
743  sq_norm_oldW = sq_norm_W;
744  start_time = get_time();
745  compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
746  ocas.w_time += get_time() - start_time;
747 
748  /* select a new cut */
749  switch( Method )
750  {
751  /* cutting plane algorithm implemented in SVMperf and BMRM */
752  case 0:
753 
754  start_time = get_time();
755  if( compute_output( output, user_data ) != 0)
756  {
757  ocas.exitflag=-2;
758  goto cleanup;
759  }
760  ocas.output_time += get_time()-start_time;
761  gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P);
762  SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6);
763 
764  xi = 0;
765  cut_length = 0;
766  ocas.trn_err = 0;
767  for(i=0; i < nData; i++)
768  {
769  if(output[i] <= 0) ocas.trn_err++;
770 
771  if(output[i] <= 1) {
772  xi += 1 - output[i];
773  new_cut[cut_length] = i;
774  cut_length++;
775  }
776  }
777  ocas.Q_P = 0.5*sq_norm_W + C*xi;
778 
779  ocas.ocas_time = get_time() - ocas_start_time;
780 
781  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
782  ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
783  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
784  */
785 
786  start_time = get_time();
787  ocas_print(ocas);
788  ocas.print_time += get_time() - start_time;
789 
790  break;
791 
792 
793  /* Ocas strategy */
794  case 1:
795 
796  /* Linesearch */
797  A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
798  B0 = dot_prod_WoldW - sq_norm_oldW;
799 
800  memcpy( old_output, output, sizeof(float64_t)*nData );
801 
802  start_time = get_time();
803  if( compute_output( output, user_data ) != 0)
804  {
805  ocas.exitflag=-2;
806  goto cleanup;
807  }
808  ocas.output_time += get_time()-start_time;
809 
810  uint32_t num_hp = 0;
811  GradVal = B0;
812  for(i=0; i< nData; i++) {
813 
814  Ci[i] = C*(1-old_output[i]);
815  Bi[i] = C*(old_output[i] - output[i]);
816 
817  float64_t val;
818  if(Bi[i] != 0)
819  val = -Ci[i]/Bi[i];
820  else
821  val = -LIBOCAS_PLUS_INF;
822 
823  if (val>0)
824  {
825 /* hpi[num_hp] = i;*/
826  hpb[num_hp] = Bi[i];
827  hpf[num_hp] = val;
828  num_hp++;
829  }
830 
831  if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
832  GradVal += Bi[i];
833 
834  }
835 
836  t = 0;
837  if( GradVal < 0 )
838  {
839  start_time = get_time();
840 /* if( sort(hpf, hpi, num_hp) != 0)*/
841  if( sort(hpf, hpb, num_hp) != 0 )
842  {
843  ocas.exitflag=-2;
844  goto cleanup;
845  }
846  ocas.sort_time += get_time() - start_time;
847 
848  float64_t t_new, GradVal_new;
849  i = 0;
850  while( GradVal < 0 && i < num_hp )
851  {
852  t_new = hpf[i];
853  GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
854 
855  if( GradVal_new >= 0 )
856  {
857  t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
858  }
859  else
860  {
861  t = t_new;
862  i++;
863  }
864 
865  GradVal = GradVal_new;
866  }
867  }
868 
869  /*
870  t = hpf[0] - 1;
871  i = 0;
872  GradVal = t*A0 + Bsum;
873  while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) {
874  t = hpf[i];
875  Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]);
876  GradVal = t*A0 + Bsum;
877  i++;
878  }
879  */
880  t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
881 
882  t1 = t; /* new (best so far) W */
883  t2 = t+MU*(1.0-t); /* new cutting plane */
884  /* t2 = t+(1.0-t)/10.0; */
885 
886  /* update W to be the best so far solution */
887  sq_norm_W = update_W( t1, user_data );
888 
889  /* select a new cut */
890  xi = 0;
891  cut_length = 0;
892  ocas.trn_err = 0;
893  for(i=0; i < nData; i++ ) {
894 
895  if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 )
896  {
897  new_cut[cut_length] = i;
898  cut_length++;
899  }
900 
901  output[i] = old_output[i]*(1-t1) + t1*output[i];
902 
903  if( output[i] <= 1) xi += 1-output[i];
904  if( output[i] <= 0) ocas.trn_err++;
905 
906  }
907 
908  ocas.Q_P = 0.5*sq_norm_W + C*xi;
909 
910  ocas.ocas_time = get_time() - ocas_start_time;
911 
912  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
913  ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
914  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
915  */
916 
917  start_time = get_time();
918  ocas_print(ocas);
919  ocas.print_time += get_time() - start_time;
920 
921  break;
922  }
923 
924  /* Stopping conditions */
925  if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
926  if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
927  if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
928  if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
929  if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
930 
931  } /* end of the main loop */
932 
933 cleanup:
934 
935  LIBOCAS_FREE(H);
936  LIBOCAS_FREE(b);
937  LIBOCAS_FREE(alpha);
938  LIBOCAS_FREE(new_cut);
939  LIBOCAS_FREE(I);
940  LIBOCAS_FREE(diag_H);
941  LIBOCAS_FREE(output);
942  LIBOCAS_FREE(old_output);
943  LIBOCAS_FREE(hpf);
944 /* LIBOCAS_FREE(hpi);*/
945  LIBOCAS_FREE(hpb);
946  LIBOCAS_FREE(Ci);
947  LIBOCAS_FREE(Bi);
948 
949  ocas.ocas_time = get_time() - ocas_start_time;
950 
951  return(ocas);
952 }
953 
954 
955 /*----------------------------------------------------------------------
956  Binary linear Ocas-SVM solver which allows using different C for each
957  training example.
958  ----------------------------------------------------------------------*/
959 ocas_return_value_T svm_ocas_solver_difC(
960  float64_t *C,
961  uint32_t nData,
962  float64_t TolRel,
963  float64_t TolAbs,
964  float64_t QPBound,
965  float64_t MaxTime,
966  uint32_t _BufSize,
967  uint8_t Method,
968  void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
969  float64_t (*update_W)(float64_t, void*),
970  int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*),
971  int (*compute_output)(float64_t*, void* ),
972  int (*sort)(float64_t*, float64_t*, uint32_t),
973  void (*ocas_print)(ocas_return_value_T),
974  void* user_data)
975 {
976  ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
977  float64_t *b, *alpha, *diag_H;
978  float64_t *output, *old_output;
979  float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
980  float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
981  float64_t start_time, ocas_start_time;
982  float64_t qp_b = 1.0;
983  float64_t new_b;
984  uint32_t cut_length;
985  uint32_t i, *new_cut;
986  uint32_t *I;
987  uint8_t S = 1;
988  libqp_state_T qp_exitflag;
989 
990  ocas_start_time = get_time();
991  ocas.qp_solver_time = 0;
992  ocas.output_time = 0;
993  ocas.sort_time = 0;
994  ocas.add_time = 0;
995  ocas.w_time = 0;
996  ocas.print_time = 0;
997 
998  BufSize = _BufSize;
999 
1000  QPSolverTolRel = TolRel*0.5;
1001 
1002  H=NULL;
1003  b=NULL;
1004  alpha=NULL;
1005  new_cut=NULL;
1006  I=NULL;
1007  diag_H=NULL;
1008  output=NULL;
1009  old_output=NULL;
1010  hpf=NULL;
1011  hpb = NULL;
1012  Ci=NULL;
1013  Bi=NULL;
1014 
1015  /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
1016  H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
1017  if(H == NULL)
1018  {
1019  ocas.exitflag=-2;
1020  goto cleanup;
1021  }
1022 
1023  /* bias of cutting planes */
1024  b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1025  if(b == NULL)
1026  {
1027  ocas.exitflag=-2;
1028  goto cleanup;
1029  }
1030 
1031  alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1032  if(alpha == NULL)
1033  {
1034  ocas.exitflag=-2;
1035  goto cleanup;
1036  }
1037 
1038  /* indices of examples which define a new cut */
1039  new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
1040  if(new_cut == NULL)
1041  {
1042  ocas.exitflag=-2;
1043  goto cleanup;
1044  }
1045 
1046  I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
1047  if(I == NULL)
1048  {
1049  ocas.exitflag=-2;
1050  goto cleanup;
1051  }
1052 
1053  for(i=0; i< BufSize; i++) I[i] = 1;
1054 
1055  diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1056  if(diag_H == NULL)
1057  {
1058  ocas.exitflag=-2;
1059  goto cleanup;
1060  }
1061 
1062  output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
1063  if(output == NULL)
1064  {
1065  ocas.exitflag=-2;
1066  goto cleanup;
1067  }
1068 
1069  old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
1070  if(old_output == NULL)
1071  {
1072  ocas.exitflag=-2;
1073  goto cleanup;
1074  }
1075 
1076  /* array of hinge points used in line-serach */
1077  hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
1078  if(hpf == NULL)
1079  {
1080  ocas.exitflag=-2;
1081  goto cleanup;
1082  }
1083 
1084  hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
1085  if(hpb == NULL)
1086  {
1087  ocas.exitflag=-2;
1088  goto cleanup;
1089  }
1090 
1091  /* vectors Ci, Bi are used in the line search procedure */
1092  Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
1093  if(Ci == NULL)
1094  {
1095  ocas.exitflag=-2;
1096  goto cleanup;
1097  }
1098 
1099  Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
1100  if(Bi == NULL)
1101  {
1102  ocas.exitflag=-2;
1103  goto cleanup;
1104  }
1105 
1106  ocas.nCutPlanes = 0;
1107  ocas.exitflag = 0;
1108  ocas.nIter = 0;
1109 
1110  /* Compute initial value of Q_P assuming that W is zero vector.*/
1111  sq_norm_W = 0;
1112  xi = nData;
1113 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/
1114  ocas.Q_D = 0;
1115 
1116  /* Compute the initial cutting plane */
1117  cut_length = nData;
1118  new_b = 0;
1119  for(i=0; i < nData; i++)
1120  {
1121  new_cut[i] = i;
1122  new_b += C[i];
1123  }
1124 
1125  ocas.Q_P = 0.5*sq_norm_W + new_b;
1126 
1127 
1128  ocas.trn_err = nData;
1129  ocas.ocas_time = get_time() - ocas_start_time;
1130  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n",
1131  ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P));
1132  */
1133  ocas_print(ocas);
1134 
1135  /* main loop */
1136  while( ocas.exitflag == 0 )
1137  {
1138  ocas.nIter++;
1139 
1140  /* append a new cut to the buffer and update H */
1141 /* b[ocas.nCutPlanes] = -(float64_t)cut_length*C;*/
1142  b[ocas.nCutPlanes] = -new_b;
1143 
1144  start_time = get_time();
1145 
1146  if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
1147  {
1148  ocas.exitflag=-2;
1149  goto cleanup;
1150  }
1151 
1152  ocas.add_time += get_time() - start_time;
1153 
1154  /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
1155  diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
1156  for(i=0; i < ocas.nCutPlanes; i++) {
1157  H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
1158  }
1159 
1160  ocas.nCutPlanes++;
1161 
1162  /* call inner QP solver */
1163  start_time = get_time();
1164 
1165 /* qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,*/
1166 /* ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);*/
1167  qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &qp_b, I, &S, alpha,
1168  ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
1169 
1170  ocas.qp_exitflag = qp_exitflag.exitflag;
1171 
1172  ocas.qp_solver_time += get_time() - start_time;
1173  ocas.Q_D = -qp_exitflag.QP;
1174 
1175  ocas.nNZAlpha = 0;
1176  for(i=0; i < ocas.nCutPlanes; i++) {
1177  if( alpha[i] != 0) ocas.nNZAlpha++;
1178  }
1179 
1180  sq_norm_oldW = sq_norm_W;
1181  start_time = get_time();
1182  compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
1183  ocas.w_time += get_time() - start_time;
1184 
1185  /* select a new cut */
1186  switch( Method )
1187  {
1188  /* cutting plane algorithm implemented in SVMperf and BMRM */
1189  case 0:
1190 
1191  start_time = get_time();
1192  if( compute_output( output, user_data ) != 0)
1193  {
1194  ocas.exitflag=-2;
1195  goto cleanup;
1196  }
1197  ocas.output_time += get_time()-start_time;
1198 
1199  xi = 0;
1200  cut_length = 0;
1201  ocas.trn_err = 0;
1202  new_b = 0;
1203  for(i=0; i < nData; i++)
1204  {
1205  if(output[i] <= 0) ocas.trn_err++;
1206 
1207 /* if(output[i] <= 1) {*/
1208 /* xi += 1 - output[i];*/
1209  if(output[i] <= C[i]) {
1210  xi += C[i] - output[i];
1211  new_cut[cut_length] = i;
1212  cut_length++;
1213  new_b += C[i];
1214  }
1215  }
1216 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/
1217  ocas.Q_P = 0.5*sq_norm_W + xi;
1218 
1219  ocas.ocas_time = get_time() - ocas_start_time;
1220 
1221  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
1222  ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
1223  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
1224  */
1225 
1226  start_time = get_time();
1227  ocas_print(ocas);
1228  ocas.print_time += get_time() - start_time;
1229 
1230  break;
1231 
1232 
1233  /* Ocas strategy */
1234  case 1:
1235 
1236  /* Linesearch */
1237  A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
1238  B0 = dot_prod_WoldW - sq_norm_oldW;
1239 
1240  memcpy( old_output, output, sizeof(float64_t)*nData );
1241 
1242  start_time = get_time();
1243  if( compute_output( output, user_data ) != 0)
1244  {
1245  ocas.exitflag=-2;
1246  goto cleanup;
1247  }
1248  ocas.output_time += get_time()-start_time;
1249 
1250  uint32_t num_hp = 0;
1251  GradVal = B0;
1252  for(i=0; i< nData; i++) {
1253 
1254 /* Ci[i] = C*(1-old_output[i]);*/
1255 /* Bi[i] = C*(old_output[i] - output[i]);*/
1256  Ci[i] = (C[i]-old_output[i]);
1257  Bi[i] = old_output[i] - output[i];
1258 
1259  float64_t val;
1260  if(Bi[i] != 0)
1261  val = -Ci[i]/Bi[i];
1262  else
1263  val = -LIBOCAS_PLUS_INF;
1264 
1265  if (val>0)
1266  {
1267 /* hpi[num_hp] = i;*/
1268  hpb[num_hp] = Bi[i];
1269  hpf[num_hp] = val;
1270  num_hp++;
1271  }
1272 
1273  if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
1274  GradVal += Bi[i];
1275 
1276  }
1277 
1278  t = 0;
1279  if( GradVal < 0 )
1280  {
1281  start_time = get_time();
1282 /* if( sort(hpf, hpi, num_hp) != 0)*/
1283  if( sort(hpf, hpb, num_hp) != 0 )
1284  {
1285  ocas.exitflag=-2;
1286  goto cleanup;
1287  }
1288  ocas.sort_time += get_time() - start_time;
1289 
1290  float64_t t_new, GradVal_new;
1291  i = 0;
1292  while( GradVal < 0 && i < num_hp )
1293  {
1294  t_new = hpf[i];
1295  GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
1296 
1297  if( GradVal_new >= 0 )
1298  {
1299  t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
1300  }
1301  else
1302  {
1303  t = t_new;
1304  i++;
1305  }
1306 
1307  GradVal = GradVal_new;
1308  }
1309  }
1310 
1311  /*
1312  t = hpf[0] - 1;
1313  i = 0;
1314  GradVal = t*A0 + Bsum;
1315  while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) {
1316  t = hpf[i];
1317  Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]);
1318  GradVal = t*A0 + Bsum;
1319  i++;
1320  }
1321  */
1322  t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
1323 
1324  t1 = t; /* new (best so far) W */
1325  t2 = t+(1.0-t)*MU; /* new cutting plane */
1326  /* t2 = t+(1.0-t)/10.0; new cutting plane */
1327 
1328  /* update W to be the best so far solution */
1329  sq_norm_W = update_W( t1, user_data );
1330 
1331  /* select a new cut */
1332  xi = 0;
1333  cut_length = 0;
1334  ocas.trn_err = 0;
1335  new_b = 0;
1336  for(i=0; i < nData; i++ ) {
1337 
1338 /* if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 ) */
1339  if( (old_output[i]*(1-t2) + t2*output[i]) <= C[i] )
1340  {
1341  new_cut[cut_length] = i;
1342  cut_length++;
1343  new_b += C[i];
1344  }
1345 
1346  output[i] = old_output[i]*(1-t1) + t1*output[i];
1347 
1348 /* if( output[i] <= 1) xi += 1-output[i];*/
1349  if( output[i] <= C[i]) xi += C[i]-output[i];
1350  if( output[i] <= 0) ocas.trn_err++;
1351 
1352  }
1353 
1354 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/
1355  ocas.Q_P = 0.5*sq_norm_W + xi;
1356 
1357  ocas.ocas_time = get_time() - ocas_start_time;
1358 
1359  /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
1360  ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
1361  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
1362  */
1363 
1364  start_time = get_time();
1365  ocas_print(ocas);
1366  ocas.print_time += get_time() - start_time;
1367 
1368  break;
1369  }
1370 
1371  /* Stopping conditions */
1372  if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
1373  if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
1374  if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
1375  if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
1376  if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
1377 
1378  } /* end of the main loop */
1379 
1380 cleanup:
1381 
1382  LIBOCAS_FREE(H);
1383  LIBOCAS_FREE(b);
1384  LIBOCAS_FREE(alpha);
1385  LIBOCAS_FREE(new_cut);
1386  LIBOCAS_FREE(I);
1387  LIBOCAS_FREE(diag_H);
1388  LIBOCAS_FREE(output);
1389  LIBOCAS_FREE(old_output);
1390  LIBOCAS_FREE(hpf);
1391 /* LIBOCAS_FREE(hpi);*/
1392  LIBOCAS_FREE(hpb);
1393  LIBOCAS_FREE(Ci);
1394  LIBOCAS_FREE(Bi);
1395 
1396  ocas.ocas_time = get_time() - ocas_start_time;
1397 
1398  return(ocas);
1399 }
1400 
1401 
1402 
1403 /*----------------------------------------------------------------------
1404  Multiclass SVM-Ocas solver
1405  ----------------------------------------------------------------------*/
1406 
1407 /* Helper function needed by the multi-class SVM linesearch.
1408 
1409  - This function finds a simplified representation of a piece-wise linear function
1410  by splitting the domain into intervals and fining active terms for these intevals */
1411 static void findactive(float64_t *Theta, float64_t *SortedA, uint32_t *nSortedA, float64_t *A, float64_t *B, int n,
1412  int (*sort)(float64_t*, float64_t*, uint32_t))
1413 {
1414  float64_t tmp, theta;
1415  uint32_t i, j, idx, idx2 = 0, start;
1416 
1417  sort(A,B,n);
1418 
1419  tmp = B[0];
1420  idx = 0;
1421  i = 0;
1422  while( i < (uint32_t)n-1 && A[i] == A[i+1])
1423  {
1424  if( B[i+1] > B[idx] )
1425  {
1426  idx = i+1;
1427  tmp = B[i+1];
1428  }
1429  i++;
1430  }
1431 
1432  (*nSortedA) = 1;
1433  SortedA[0] = A[idx];
1434 
1435  while(1)
1436  {
1437  start = idx + 1;
1438  while( start < (uint32_t)n && A[idx] == A[start])
1439  start++;
1440 
1441  theta = LIBOCAS_PLUS_INF;
1442  for(j=start; j < (uint32_t)n; j++)
1443  {
1444  tmp = (B[j] - B[idx])/(A[idx]-A[j]);
1445  if( tmp < theta)
1446  {
1447  theta = tmp;
1448  idx2 = j;
1449  }
1450  }
1451 
1452  if( theta < LIBOCAS_PLUS_INF)
1453  {
1454  Theta[(*nSortedA) - 1] = theta;
1455  SortedA[(*nSortedA)] = A[idx2];
1456  (*nSortedA)++;
1457  idx = idx2;
1458  }
1459  else
1460  return;
1461  }
1462 }
1463 
1464 
1465 /*----------------------------------------------------------------------
1466  Multiclass linear OCAS-SVM solver.
1467  ----------------------------------------------------------------------*/
1468 ocas_return_value_T msvm_ocas_solver(
1469  float64_t C,
1470  float64_t *data_y,
1471  uint32_t nY,
1472  uint32_t nData,
1473  float64_t TolRel,
1474  float64_t TolAbs,
1475  float64_t QPBound,
1476  float64_t MaxTime,
1477  uint32_t _BufSize,
1478  uint8_t Method,
1479  void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
1480  float64_t (*update_W)(float64_t, void*),
1481  int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, void*),
1482  int (*compute_output)(float64_t*, void* ),
1483  int (*sort)(float64_t*, float64_t*, uint32_t),
1484  void (*ocas_print)(ocas_return_value_T),
1485  void* user_data)
1486 {
1487  ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1488  float64_t *b, *alpha, *diag_H;
1489  float64_t *output, *old_output;
1490  float64_t xi, sq_norm_W, QPSolverTolRel, QPSolverTolAbs, dot_prod_WoldW, sq_norm_oldW;
1491  float64_t A0, B0, t, t1, t2, R, tmp, element_b, x;
1492  float64_t *A, *B, *theta, *Theta, *sortedA, *Add;
1493  float64_t start_time, ocas_start_time, grad_sum, grad, min_x = 0, old_x, old_grad;
1494  uint32_t i, y, y2, ypred = 0, *new_cut, cnt1, cnt2, j, nSortedA, idx;
1495  uint32_t *I;
1496  uint8_t S = 1;
1497  libqp_state_T qp_exitflag;
1498 
1499  ocas_start_time = get_time();
1500  ocas.qp_solver_time = 0;
1501  ocas.output_time = 0;
1502  ocas.sort_time = 0;
1503  ocas.add_time = 0;
1504  ocas.w_time = 0;
1505  ocas.print_time = 0;
1506 
1507  BufSize = _BufSize;
1508 
1509  QPSolverTolRel = TolRel*0.5;
1510  QPSolverTolAbs = TolAbs*0.5;
1511 
1512  H=NULL;
1513  b=NULL;
1514  alpha=NULL;
1515  new_cut=NULL;
1516  I=NULL;
1517  diag_H=NULL;
1518  output=NULL;
1519  old_output=NULL;
1520  A = NULL;
1521  B = NULL;
1522  theta = NULL;
1523  Theta = NULL;
1524  sortedA = NULL;
1525  Add = NULL;
1526 
1527  /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
1528  H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
1529  if(H == NULL)
1530  {
1531  ocas.exitflag=-2;
1532  goto cleanup;
1533  }
1534 
1535  /* bias of cutting planes */
1536  b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1537  if(b == NULL)
1538  {
1539  ocas.exitflag=-2;
1540  goto cleanup;
1541  }
1542 
1543  alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1544  if(alpha == NULL)
1545  {
1546  ocas.exitflag=-2;
1547  goto cleanup;
1548  }
1549 
1550  /* indices of examples which define a new cut */
1551  new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
1552  if(new_cut == NULL)
1553  {
1554  ocas.exitflag=-2;
1555  goto cleanup;
1556  }
1557 
1558  I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
1559  if(I == NULL)
1560  {
1561  ocas.exitflag=-2;
1562  goto cleanup;
1563  }
1564 
1565  for(i=0; i< BufSize; i++)
1566  I[i] = 1;
1567 
1568  diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1569  if(diag_H == NULL)
1570  {
1571  ocas.exitflag=-2;
1572  goto cleanup;
1573  }
1574 
1575  output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1576  if(output == NULL)
1577  {
1578  ocas.exitflag=-2;
1579  goto cleanup;
1580  }
1581 
1582  old_output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1583  if(old_output == NULL)
1584  {
1585  ocas.exitflag=-2;
1586  goto cleanup;
1587  }
1588 
1589  /* auxciliary variables used in the linesearch */
1590  A = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1591  if(A == NULL)
1592  {
1593  ocas.exitflag=-2;
1594  goto cleanup;
1595  }
1596 
1597  B = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1598  if(B == NULL)
1599  {
1600  ocas.exitflag=-2;
1601  goto cleanup;
1602  }
1603 
1604  theta = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t));
1605  if(theta == NULL)
1606  {
1607  ocas.exitflag=-2;
1608  goto cleanup;
1609  }
1610 
1611  sortedA = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t));
1612  if(sortedA == NULL)
1613  {
1614  ocas.exitflag=-2;
1615  goto cleanup;
1616  }
1617 
1618  Theta = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1619  if(Theta == NULL)
1620  {
1621  ocas.exitflag=-2;
1622  goto cleanup;
1623  }
1624 
1625  Add = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1626  if(Add == NULL)
1627  {
1628  ocas.exitflag=-2;
1629  goto cleanup;
1630  }
1631 
1632  /* Set initial values*/
1633  ocas.nCutPlanes = 0;
1634  ocas.exitflag = 0;
1635  ocas.nIter = 0;
1636  ocas.Q_D = 0;
1637  ocas.trn_err = nData;
1638  R = (float64_t)nData;
1639  sq_norm_W = 0;
1640  element_b = (float64_t)nData;
1641  ocas.Q_P = 0.5*sq_norm_W + C*R;
1642 
1643  /* initial cutting plane */
1644  for(i=0; i < nData; i++)
1645  {
1646  y2 = (uint32_t)data_y[i];
1647 
1648  if(y2 > 0)
1649  new_cut[i] = 0;
1650  else
1651  new_cut[i] = 1;
1652 
1653  }
1654 
1655  ocas.ocas_time = get_time() - ocas_start_time;
1656 
1657  start_time = get_time();
1658  ocas_print(ocas);
1659  ocas.print_time += get_time() - start_time;
1660 
1661  /* main loop of the OCAS */
1662  while( ocas.exitflag == 0 )
1663  {
1664  ocas.nIter++;
1665 
1666  /* append a new cut to the buffer and update H */
1667  b[ocas.nCutPlanes] = -(float64_t)element_b;
1668 
1669  start_time = get_time();
1670 
1671  if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, ocas.nCutPlanes, user_data ) != 0)
1672  {
1673  ocas.exitflag=-2;
1674  goto cleanup;
1675  }
1676 
1677  ocas.add_time += get_time() - start_time;
1678 
1679  /* copy newly appended row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
1680  diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
1681  for(i=0; i < ocas.nCutPlanes; i++)
1682  {
1683  H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
1684  }
1685 
1686  ocas.nCutPlanes++;
1687 
1688  /* call inner QP solver */
1689  start_time = get_time();
1690 
1691  qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
1692  ocas.nCutPlanes, QPSolverMaxIter, QPSolverTolAbs, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
1693 
1694  ocas.qp_exitflag = qp_exitflag.exitflag;
1695 
1696  ocas.qp_solver_time += get_time() - start_time;
1697  ocas.Q_D = -qp_exitflag.QP;
1698 
1699  ocas.nNZAlpha = 0;
1700  for(i=0; i < ocas.nCutPlanes; i++)
1701  if( alpha[i] != 0) ocas.nNZAlpha++;
1702 
1703  sq_norm_oldW = sq_norm_W;
1704  start_time = get_time();
1705  compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
1706  ocas.w_time += get_time() - start_time;
1707 
1708  /* select a new cut */
1709  switch( Method )
1710  {
1711  /* cutting plane algorithm implemented in SVMperf and BMRM */
1712  case 0:
1713 
1714  start_time = get_time();
1715  if( compute_output( output, user_data ) != 0)
1716  {
1717  ocas.exitflag=-2;
1718  goto cleanup;
1719  }
1720  ocas.output_time += get_time()-start_time;
1721 
1722  /* the following loop computes: */
1723  element_b = 0.0; /* element_b = R(old_W) - g'*old_W */
1724  R = 0; /* R(W) = sum_i max_y ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */
1725  ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */
1726  /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */
1727  for(i=0; i < nData; i++)
1728  {
1729  y2 = (uint32_t)data_y[i];
1730 
1731  for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1732  {
1733  if(y2 != y && xi < output[LIBOCAS_INDEX(y,i,nY)])
1734  {
1735  xi = output[LIBOCAS_INDEX(y,i,nY)];
1736  ypred = y;
1737  }
1738  }
1739 
1740  if(xi >= output[LIBOCAS_INDEX(y2,i,nY)])
1741  ocas.trn_err ++;
1742 
1743  xi = LIBOCAS_MAX(0,xi+1-output[LIBOCAS_INDEX(y2,i,nY)]);
1744  R += xi;
1745  if(xi > 0)
1746  {
1747  element_b++;
1748  new_cut[i] = ypred;
1749  }
1750  else
1751  new_cut[i] = y2;
1752  }
1753 
1754  ocas.Q_P = 0.5*sq_norm_W + C*R;
1755 
1756  ocas.ocas_time = get_time() - ocas_start_time;
1757 
1758  start_time = get_time();
1759  ocas_print(ocas);
1760  ocas.print_time += get_time() - start_time;
1761 
1762  break;
1763 
1764  /* The OCAS solver */
1765  case 1:
1766  memcpy( old_output, output, sizeof(float64_t)*nData*nY );
1767 
1768  start_time = get_time();
1769  if( compute_output( output, user_data ) != 0)
1770  {
1771  ocas.exitflag=-2;
1772  goto cleanup;
1773  }
1774  ocas.output_time += get_time()-start_time;
1775 
1776  A0 = sq_norm_W - 2*dot_prod_WoldW + sq_norm_oldW;
1777  B0 = dot_prod_WoldW - sq_norm_oldW;
1778 
1779  for(i=0; i < nData; i++)
1780  {
1781  y2 = (uint32_t)data_y[i];
1782 
1783  for(y=0; y < nY; y++)
1784  {
1785  A[LIBOCAS_INDEX(y,i,nY)] = C*(output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y,i,nY)]
1786  + old_output[LIBOCAS_INDEX(y2,i,nY)] - output[LIBOCAS_INDEX(y2,i,nY)]);
1787  B[LIBOCAS_INDEX(y,i,nY)] = C*(old_output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y2,i,nY)]
1788  + (float64_t)(y != y2));
1789  }
1790  }
1791 
1792  /* linesearch */
1793 /* new_x = msvm_linesearch_mex(A0,B0,AA*C,BB*C);*/
1794 
1795  grad_sum = B0;
1796  cnt1 = 0;
1797  cnt2 = 0;
1798  for(i=0; i < nData; i++)
1799  {
1800  findactive(theta,sortedA,&nSortedA,&A[i*nY],&B[i*nY],nY,sort);
1801 
1802  idx = 0;
1803  while( idx < nSortedA-1 && theta[idx] < 0 )
1804  idx++;
1805 
1806  grad_sum += sortedA[idx];
1807 
1808  for(j=idx; j < nSortedA-1; j++)
1809  {
1810  Theta[cnt1] = theta[j];
1811  cnt1++;
1812  }
1813 
1814  for(j=idx+1; j < nSortedA; j++)
1815  {
1816  Add[cnt2] = -sortedA[j-1]+sortedA[j];
1817  cnt2++;
1818  }
1819  }
1820 
1821  start_time = get_time();
1822  sort(Theta,Add,cnt1);
1823  ocas.sort_time += get_time() - start_time;
1824 
1825  grad = grad_sum;
1826  if(grad >= 0)
1827  {
1828  min_x = 0;
1829  }
1830  else
1831  {
1832  old_x = 0;
1833  old_grad = grad;
1834 
1835  for(i=0; i < cnt1; i++)
1836  {
1837  x = Theta[i];
1838 
1839  grad = x*A0 + grad_sum;
1840 
1841  if(grad >=0)
1842  {
1843 
1844  min_x = (grad*old_x - old_grad*x)/(grad - old_grad);
1845 
1846  break;
1847  }
1848  else
1849  {
1850  grad_sum = grad_sum + Add[i];
1851 
1852  grad = x*A0 + grad_sum;
1853  if( grad >= 0)
1854  {
1855  min_x = x;
1856  break;
1857  }
1858  }
1859 
1860  old_grad = grad;
1861  old_x = x;
1862  }
1863  }
1864  /* end of the linesearch which outputs min_x */
1865 
1866  t = min_x;
1867  t1 = t; /* new (best so far) W */
1868  t2 = t+(1.0-t)*MU; /* new cutting plane */
1869  /* t2 = t+(1.0-t)/10.0; */
1870 
1871  /* update W to be the best so far solution */
1872  sq_norm_W = update_W( t1, user_data );
1873 
1874  /* the following code computes a new cutting plane: */
1875  element_b = 0.0; /* element_b = R(old_W) - g'*old_W */
1876  /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */
1877  for(i=0; i < nData; i++)
1878  {
1879  y2 = (uint32_t)data_y[i];
1880 
1881  for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1882  {
1883  tmp = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y,i,nY)];
1884  if(y2 != y && xi < tmp)
1885  {
1886  xi = tmp;
1887  ypred = y;
1888  }
1889  }
1890 
1891  tmp = old_output[LIBOCAS_INDEX(y2,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y2,i,nY)];
1892  xi = LIBOCAS_MAX(0,xi+1-tmp);
1893  if(xi > 0)
1894  {
1895  element_b++;
1896  new_cut[i] = ypred;
1897  }
1898  else
1899  new_cut[i] = y2;
1900  }
1901 
1902  /* compute Risk, class. error and update outputs to correspond to the new W */
1903  ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */
1904  R = 0;
1905  for(i=0; i < nData; i++)
1906  {
1907  y2 = (uint32_t)data_y[i];
1908 
1909  for(tmp=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1910  {
1911  output[LIBOCAS_INDEX(y,i,nY)] = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t1) + t1*output[LIBOCAS_INDEX(y,i,nY)];
1912 
1913  if(y2 != y && tmp < output[LIBOCAS_INDEX(y,i,nY)])
1914  {
1915  ypred = y;
1916  tmp = output[LIBOCAS_INDEX(y,i,nY)];
1917  }
1918  }
1919 
1920  R += LIBOCAS_MAX(0,1+tmp - output[LIBOCAS_INDEX(y2,i,nY)]);
1921  if( tmp >= output[LIBOCAS_INDEX(y2,i,nY)])
1922  ocas.trn_err ++;
1923  }
1924 
1925  ocas.Q_P = 0.5*sq_norm_W + C*R;
1926 
1927 
1928  /* get time and print status */
1929  ocas.ocas_time = get_time() - ocas_start_time;
1930 
1931  start_time = get_time();
1932  ocas_print(ocas);
1933  ocas.print_time += get_time() - start_time;
1934 
1935  break;
1936 
1937  }
1938 
1939  /* Stopping conditions */
1940  if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
1941  if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
1942  if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
1943  if( MaxTime > 0 && ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
1944  if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
1945 
1946  } /* end of the main loop */
1947 
1948 cleanup:
1949 
1950  LIBOCAS_FREE(H);
1951  LIBOCAS_FREE(b);
1952  LIBOCAS_FREE(alpha);
1953  LIBOCAS_FREE(new_cut);
1954  LIBOCAS_FREE(I);
1955  LIBOCAS_FREE(diag_H);
1956  LIBOCAS_FREE(output);
1957  LIBOCAS_FREE(old_output);
1958  LIBOCAS_FREE(A);
1959  LIBOCAS_FREE(B);
1960  LIBOCAS_FREE(theta);
1961  LIBOCAS_FREE(Theta);
1962  LIBOCAS_FREE(sortedA);
1963  LIBOCAS_FREE(Add);
1964 
1965  ocas.ocas_time = get_time() - ocas_start_time;
1966 
1967  return(ocas);
1968 }
1969 }
1970 
1971 

SHOGUN Machine Learning Toolbox - Documentation