SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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 LAMBDA 0.1 /* must be from (0,1> 1..means that OCAS becomes equivalent to CPA */
29 
30 static const uint32_t QPSolverMaxIter = 10000000;
31 
32 static float64_t *H;
33 static uint32_t BufSize;
34 
35 /*----------------------------------------------------------------------
36  Returns pointer at i-th column of Hessian matrix.
37  ----------------------------------------------------------------------*/
38 static const float64_t *get_col( uint32_t i)
39 {
40  return( &H[ BufSize*i ] );
41 }
42 
43 /*----------------------------------------------------------------------
44  Returns time of the day in seconds.
45  ----------------------------------------------------------------------*/
47 {
48  struct timeval tv;
49  if (gettimeofday(&tv, NULL)==0)
50  return tv.tv_sec+((float64_t)(tv.tv_usec))/1e6;
51  else
52  return 0.0;
53 }
54 
55 /*----------------------------------------------------------------------
56  Linear binary Ocas-SVM solver.
57  ----------------------------------------------------------------------*/
58 ocas_return_value_T svm_ocas_solver(
59  float64_t C,
60  uint32_t nData,
61  float64_t TolRel,
62  float64_t TolAbs,
63  float64_t QPBound,
64  float64_t MaxTime,
65  uint32_t _BufSize,
66  uint8_t Method,
67  void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
68  float64_t (*update_W)(float64_t, void*),
69  int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*),
70  int (*compute_output)(float64_t*, void* ),
71  int (*sort)(float64_t*, float64_t*, uint32_t),
72  void (*ocas_print)(ocas_return_value_T),
73  void* user_data)
74 {
75  ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
76  float64_t *b, *alpha, *diag_H;
77  float64_t *output, *old_output;
78  float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
79  float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
80  float64_t start_time, ocas_start_time;
81  uint32_t cut_length;
82  uint32_t i, *new_cut;
83  uint32_t *I;
84  uint8_t S = 1;
85  libqp_state_T qp_exitflag;
86 
87  ocas_start_time = get_time();
88  ocas.qp_solver_time = 0;
89  ocas.output_time = 0;
90  ocas.sort_time = 0;
91  ocas.add_time = 0;
92  ocas.w_time = 0;
93  ocas.print_time = 0;
94  float64_t gap;
95 
96  BufSize = _BufSize;
97 
98  QPSolverTolRel = TolRel*0.5;
99 
100  H=NULL;
101  b=NULL;
102  alpha=NULL;
103  new_cut=NULL;
104  I=NULL;
105  diag_H=NULL;
106  output=NULL;
107  old_output=NULL;
108  hpf=NULL;
109  hpb = NULL;
110  Ci=NULL;
111  Bi=NULL;
112 
113  /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
114  H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
115  if(H == NULL)
116  {
117  ocas.exitflag=-2;
118  goto cleanup;
119  }
120 
121  /* bias of cutting planes */
122  b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
123  if(b == NULL)
124  {
125  ocas.exitflag=-2;
126  goto cleanup;
127  }
128 
129  alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
130  if(alpha == NULL)
131  {
132  ocas.exitflag=-2;
133  goto cleanup;
134  }
135 
136  /* indices of examples which define a new cut */
137  new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
138  if(new_cut == NULL)
139  {
140  ocas.exitflag=-2;
141  goto cleanup;
142  }
143 
144  I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
145  if(I == NULL)
146  {
147  ocas.exitflag=-2;
148  goto cleanup;
149  }
150 
151  for(i=0; i< BufSize; i++) I[i] = 1;
152 
153  diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
154  if(diag_H == NULL)
155  {
156  ocas.exitflag=-2;
157  goto cleanup;
158  }
159 
160  output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
161  if(output == NULL)
162  {
163  ocas.exitflag=-2;
164  goto cleanup;
165  }
166 
167  old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
168  if(old_output == NULL)
169  {
170  ocas.exitflag=-2;
171  goto cleanup;
172  }
173 
174  /* array of hinge points used in line-serach */
175  hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
176  if(hpf == NULL)
177  {
178  ocas.exitflag=-2;
179  goto cleanup;
180  }
181 
182  hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
183  if(hpb == NULL)
184  {
185  ocas.exitflag=-2;
186  goto cleanup;
187  }
188 
189  /* vectors Ci, Bi are used in the line search procedure */
190  Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
191  if(Ci == NULL)
192  {
193  ocas.exitflag=-2;
194  goto cleanup;
195  }
196 
197  Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
198  if(Bi == NULL)
199  {
200  ocas.exitflag=-2;
201  goto cleanup;
202  }
203 
204  ocas.nCutPlanes = 0;
205  ocas.exitflag = 0;
206  ocas.nIter = 0;
207 
208  /* Compute initial value of Q_P assuming that W is zero vector.*/
209  sq_norm_W = 0;
210  xi = nData;
211  ocas.Q_P = 0.5*sq_norm_W + C*xi;
212  ocas.Q_D = 0;
213 
214  /* Compute the initial cutting plane */
215  cut_length = nData;
216  for(i=0; i < nData; i++)
217  new_cut[i] = i;
218 
219  gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P);
220  SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6);
221 
222  ocas.trn_err = nData;
223  ocas.ocas_time = get_time() - ocas_start_time;
224  /* 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",
225  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));
226  */
227  ocas_print(ocas);
228 
229  /* main loop */
230  while( ocas.exitflag == 0 )
231  {
232  ocas.nIter++;
233 
234  /* append a new cut to the buffer and update H */
235  b[ocas.nCutPlanes] = -(float64_t)cut_length;
236 
237  start_time = get_time();
238 
239  if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
240  {
241  ocas.exitflag=-2;
242  goto cleanup;
243  }
244 
245  ocas.add_time += get_time() - start_time;
246 
247  /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
248  diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
249  for(i=0; i < ocas.nCutPlanes; i++) {
250  H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
251  }
252 
253  ocas.nCutPlanes++;
254 
255  /* call inner QP solver */
256  start_time = get_time();
257 
258  qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
259  ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
260 
261  ocas.qp_exitflag = qp_exitflag.exitflag;
262 
263  ocas.qp_solver_time += get_time() - start_time;
264  ocas.Q_D = -qp_exitflag.QP;
265 
266  ocas.nNZAlpha = 0;
267  for(i=0; i < ocas.nCutPlanes; i++) {
268  if( alpha[i] != 0) ocas.nNZAlpha++;
269  }
270 
271  sq_norm_oldW = sq_norm_W;
272  start_time = get_time();
273  compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
274  ocas.w_time += get_time() - start_time;
275 
276  /* select a new cut */
277  switch( Method )
278  {
279  /* cutting plane algorithm implemented in SVMperf and BMRM */
280  case 0:
281 
282  start_time = get_time();
283  if( compute_output( output, user_data ) != 0)
284  {
285  ocas.exitflag=-2;
286  goto cleanup;
287  }
288  ocas.output_time += get_time()-start_time;
289 
290  xi = 0;
291  cut_length = 0;
292  ocas.trn_err = 0;
293  for(i=0; i < nData; i++)
294  {
295  if(output[i] <= 0) ocas.trn_err++;
296 
297  if(output[i] <= 1) {
298  xi += 1 - output[i];
299  new_cut[cut_length] = i;
300  cut_length++;
301  }
302  }
303  ocas.Q_P = 0.5*sq_norm_W + C*xi;
304 
305  gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P);
306  SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6);
307  /* 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",
308  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),
309  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
310  */
311 
312  start_time = get_time();
313  ocas_print(ocas);
314  ocas.print_time += get_time() - start_time;
315 
316  break;
317 
318 
319  /* Ocas strategy */
320  case 1:
321 
322  /* Linesearch */
323  A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
324  B0 = dot_prod_WoldW - sq_norm_oldW;
325 
326  memcpy( old_output, output, sizeof(float64_t)*nData );
327 
328  start_time = get_time();
329  if( compute_output( output, user_data ) != 0)
330  {
331  ocas.exitflag=-2;
332  goto cleanup;
333  }
334  ocas.output_time += get_time()-start_time;
335 
336  uint32_t num_hp = 0;
337  GradVal = B0;
338  for(i=0; i< nData; i++) {
339 
340  Ci[i] = C*(1-old_output[i]);
341  Bi[i] = C*(old_output[i] - output[i]);
342 
343  float64_t val;
344  if(Bi[i] != 0)
345  val = -Ci[i]/Bi[i];
346  else
347  val = -LIBOCAS_PLUS_INF;
348 
349  if (val>0)
350  {
351 /* hpi[num_hp] = i;*/
352  hpb[num_hp] = Bi[i];
353  hpf[num_hp] = val;
354  num_hp++;
355  }
356 
357  if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
358  GradVal += Bi[i];
359 
360  }
361 
362  t = 0;
363  if( GradVal < 0 )
364  {
365  start_time = get_time();
366 /* if( sort(hpf, hpi, num_hp) != 0)*/
367  if( sort(hpf, hpb, num_hp) != 0 )
368  {
369  ocas.exitflag=-2;
370  goto cleanup;
371  }
372  ocas.sort_time += get_time() - start_time;
373 
374  float64_t t_new, GradVal_new;
375  i = 0;
376  while( GradVal < 0 && i < num_hp )
377  {
378  t_new = hpf[i];
379  GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
380 
381  if( GradVal_new >= 0 )
382  {
383  t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
384  }
385  else
386  {
387  t = t_new;
388  i++;
389  }
390 
391  GradVal = GradVal_new;
392  }
393  }
394 
395  /*
396  t = hpf[0] - 1;
397  i = 0;
398  GradVal = t*A0 + Bsum;
399  while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) {
400  t = hpf[i];
401  Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]);
402  GradVal = t*A0 + Bsum;
403  i++;
404  }
405  */
406  t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
407 
408  t1 = t; /* new (best so far) W */
409  t2 = t+LAMBDA*(1.0-t); /* new cutting plane */
410  /* t2 = t+(1.0-t)/10.0; */
411 
412  /* update W to be the best so far solution */
413  sq_norm_W = update_W( t1, user_data );
414 
415  /* select a new cut */
416  xi = 0;
417  cut_length = 0;
418  ocas.trn_err = 0;
419  for(i=0; i < nData; i++ ) {
420 
421  if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 )
422  {
423  new_cut[cut_length] = i;
424  cut_length++;
425  }
426 
427  output[i] = old_output[i]*(1-t1) + t1*output[i];
428 
429  if( output[i] <= 1) xi += 1-output[i];
430  if( output[i] <= 0) ocas.trn_err++;
431 
432  }
433 
434  ocas.Q_P = 0.5*sq_norm_W + C*xi;
435 
436  ocas.ocas_time = get_time() - ocas_start_time;
437 
438  /* 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",
439  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),
440  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
441  */
442 
443  start_time = get_time();
444  ocas_print(ocas);
445  ocas.print_time += get_time() - start_time;
446 
447  break;
448  }
449 
450  /* Stopping conditions */
451  if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
452  if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
453  if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
454  if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
455  if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
456 
457  } /* end of the main loop */
458 
459 cleanup:
460 
461  LIBOCAS_FREE(H);
462  LIBOCAS_FREE(b);
463  LIBOCAS_FREE(alpha);
464  LIBOCAS_FREE(new_cut);
465  LIBOCAS_FREE(I);
466  LIBOCAS_FREE(diag_H);
467  LIBOCAS_FREE(output);
468  LIBOCAS_FREE(old_output);
469  LIBOCAS_FREE(hpf);
470 /* LIBOCAS_FREE(hpi);*/
471  LIBOCAS_FREE(hpb);
472  LIBOCAS_FREE(Ci);
473  LIBOCAS_FREE(Bi);
474 
475  ocas.ocas_time = get_time() - ocas_start_time;
476 
477  return(ocas);
478 }
479 
480 
481 /*----------------------------------------------------------------------
482  Binary linear Ocas-SVM solver which allows using different C for each
483  training example.
484  ----------------------------------------------------------------------*/
485 ocas_return_value_T svm_ocas_solver_difC(
486  float64_t *C,
487  uint32_t nData,
488  float64_t TolRel,
489  float64_t TolAbs,
490  float64_t QPBound,
491  float64_t MaxTime,
492  uint32_t _BufSize,
493  uint8_t Method,
494  void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
495  float64_t (*update_W)(float64_t, void*),
496  int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*),
497  int (*compute_output)(float64_t*, void* ),
498  int (*sort)(float64_t*, float64_t*, uint32_t),
499  void (*ocas_print)(ocas_return_value_T),
500  void* user_data)
501 {
502  ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
503  float64_t *b, *alpha, *diag_H;
504  float64_t *output, *old_output;
505  float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
506  float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
507  float64_t start_time, ocas_start_time;
508  float64_t qp_b = 1.0;
509  float64_t new_b;
510  uint32_t cut_length;
511  uint32_t i, *new_cut;
512  uint32_t *I;
513  uint8_t S = 1;
514  libqp_state_T qp_exitflag;
515 
516  ocas_start_time = get_time();
517  ocas.qp_solver_time = 0;
518  ocas.output_time = 0;
519  ocas.sort_time = 0;
520  ocas.add_time = 0;
521  ocas.w_time = 0;
522  ocas.print_time = 0;
523 
524  BufSize = _BufSize;
525 
526  QPSolverTolRel = TolRel*0.5;
527 
528  H=NULL;
529  b=NULL;
530  alpha=NULL;
531  new_cut=NULL;
532  I=NULL;
533  diag_H=NULL;
534  output=NULL;
535  old_output=NULL;
536  hpf=NULL;
537  hpb = NULL;
538  Ci=NULL;
539  Bi=NULL;
540 
541  /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
542  H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
543  if(H == NULL)
544  {
545  ocas.exitflag=-2;
546  goto cleanup;
547  }
548 
549  /* bias of cutting planes */
550  b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
551  if(b == NULL)
552  {
553  ocas.exitflag=-2;
554  goto cleanup;
555  }
556 
557  alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
558  if(alpha == NULL)
559  {
560  ocas.exitflag=-2;
561  goto cleanup;
562  }
563 
564  /* indices of examples which define a new cut */
565  new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
566  if(new_cut == NULL)
567  {
568  ocas.exitflag=-2;
569  goto cleanup;
570  }
571 
572  I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
573  if(I == NULL)
574  {
575  ocas.exitflag=-2;
576  goto cleanup;
577  }
578 
579  for(i=0; i< BufSize; i++) I[i] = 1;
580 
581  diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
582  if(diag_H == NULL)
583  {
584  ocas.exitflag=-2;
585  goto cleanup;
586  }
587 
588  output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
589  if(output == NULL)
590  {
591  ocas.exitflag=-2;
592  goto cleanup;
593  }
594 
595  old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
596  if(old_output == NULL)
597  {
598  ocas.exitflag=-2;
599  goto cleanup;
600  }
601 
602  /* array of hinge points used in line-serach */
603  hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
604  if(hpf == NULL)
605  {
606  ocas.exitflag=-2;
607  goto cleanup;
608  }
609 
610  hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
611  if(hpb == NULL)
612  {
613  ocas.exitflag=-2;
614  goto cleanup;
615  }
616 
617  /* vectors Ci, Bi are used in the line search procedure */
618  Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
619  if(Ci == NULL)
620  {
621  ocas.exitflag=-2;
622  goto cleanup;
623  }
624 
625  Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t));
626  if(Bi == NULL)
627  {
628  ocas.exitflag=-2;
629  goto cleanup;
630  }
631 
632  ocas.nCutPlanes = 0;
633  ocas.exitflag = 0;
634  ocas.nIter = 0;
635 
636  /* Compute initial value of Q_P assuming that W is zero vector.*/
637  sq_norm_W = 0;
638  xi = nData;
639 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/
640  ocas.Q_D = 0;
641 
642  /* Compute the initial cutting plane */
643  cut_length = nData;
644  new_b = 0;
645  for(i=0; i < nData; i++)
646  {
647  new_cut[i] = i;
648  new_b += C[i];
649  }
650 
651  ocas.Q_P = 0.5*sq_norm_W + new_b;
652 
653 
654  ocas.trn_err = nData;
655  ocas.ocas_time = get_time() - ocas_start_time;
656  /* 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",
657  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));
658  */
659  ocas_print(ocas);
660 
661  /* main loop */
662  while( ocas.exitflag == 0 )
663  {
664  ocas.nIter++;
665 
666  /* append a new cut to the buffer and update H */
667 /* b[ocas.nCutPlanes] = -(float64_t)cut_length*C;*/
668  b[ocas.nCutPlanes] = -new_b;
669 
670  start_time = get_time();
671 
672  if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
673  {
674  ocas.exitflag=-2;
675  goto cleanup;
676  }
677 
678  ocas.add_time += get_time() - start_time;
679 
680  /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
681  diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
682  for(i=0; i < ocas.nCutPlanes; i++) {
683  H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
684  }
685 
686  ocas.nCutPlanes++;
687 
688  /* call inner QP solver */
689  start_time = get_time();
690 
691 /* qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,*/
692 /* ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);*/
693  qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &qp_b, I, &S, alpha,
694  ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
695 
696  ocas.qp_exitflag = qp_exitflag.exitflag;
697 
698  ocas.qp_solver_time += get_time() - start_time;
699  ocas.Q_D = -qp_exitflag.QP;
700 
701  ocas.nNZAlpha = 0;
702  for(i=0; i < ocas.nCutPlanes; i++) {
703  if( alpha[i] != 0) ocas.nNZAlpha++;
704  }
705 
706  sq_norm_oldW = sq_norm_W;
707  start_time = get_time();
708  compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
709  ocas.w_time += get_time() - start_time;
710 
711  /* select a new cut */
712  switch( Method )
713  {
714  /* cutting plane algorithm implemented in SVMperf and BMRM */
715  case 0:
716 
717  start_time = get_time();
718  if( compute_output( output, user_data ) != 0)
719  {
720  ocas.exitflag=-2;
721  goto cleanup;
722  }
723  ocas.output_time += get_time()-start_time;
724 
725  xi = 0;
726  cut_length = 0;
727  ocas.trn_err = 0;
728  new_b = 0;
729  for(i=0; i < nData; i++)
730  {
731  if(output[i] <= 0) ocas.trn_err++;
732 
733 /* if(output[i] <= 1) {*/
734 /* xi += 1 - output[i];*/
735  if(output[i] <= C[i]) {
736  xi += C[i] - output[i];
737  new_cut[cut_length] = i;
738  cut_length++;
739  new_b += C[i];
740  }
741  }
742 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/
743  ocas.Q_P = 0.5*sq_norm_W + xi;
744 
745  ocas.ocas_time = get_time() - ocas_start_time;
746 
747  /* 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",
748  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),
749  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
750  */
751 
752  start_time = get_time();
753  ocas_print(ocas);
754  ocas.print_time += get_time() - start_time;
755 
756  break;
757 
758 
759  /* Ocas strategy */
760  case 1:
761 
762  /* Linesearch */
763  A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
764  B0 = dot_prod_WoldW - sq_norm_oldW;
765 
766  memcpy( old_output, output, sizeof(float64_t)*nData );
767 
768  start_time = get_time();
769  if( compute_output( output, user_data ) != 0)
770  {
771  ocas.exitflag=-2;
772  goto cleanup;
773  }
774  ocas.output_time += get_time()-start_time;
775 
776  uint32_t num_hp = 0;
777  GradVal = B0;
778  for(i=0; i< nData; i++) {
779 
780 /* Ci[i] = C*(1-old_output[i]);*/
781 /* Bi[i] = C*(old_output[i] - output[i]);*/
782  Ci[i] = (C[i]-old_output[i]);
783  Bi[i] = old_output[i] - output[i];
784 
785  float64_t val;
786  if(Bi[i] != 0)
787  val = -Ci[i]/Bi[i];
788  else
789  val = -LIBOCAS_PLUS_INF;
790 
791  if (val>0)
792  {
793 /* hpi[num_hp] = i;*/
794  hpb[num_hp] = Bi[i];
795  hpf[num_hp] = val;
796  num_hp++;
797  }
798 
799  if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
800  GradVal += Bi[i];
801 
802  }
803 
804  t = 0;
805  if( GradVal < 0 )
806  {
807  start_time = get_time();
808 /* if( sort(hpf, hpi, num_hp) != 0)*/
809  if( sort(hpf, hpb, num_hp) != 0 )
810  {
811  ocas.exitflag=-2;
812  goto cleanup;
813  }
814  ocas.sort_time += get_time() - start_time;
815 
816  float64_t t_new, GradVal_new;
817  i = 0;
818  while( GradVal < 0 && i < num_hp )
819  {
820  t_new = hpf[i];
821  GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
822 
823  if( GradVal_new >= 0 )
824  {
825  t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
826  }
827  else
828  {
829  t = t_new;
830  i++;
831  }
832 
833  GradVal = GradVal_new;
834  }
835  }
836 
837  /*
838  t = hpf[0] - 1;
839  i = 0;
840  GradVal = t*A0 + Bsum;
841  while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) {
842  t = hpf[i];
843  Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]);
844  GradVal = t*A0 + Bsum;
845  i++;
846  }
847  */
848  t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
849 
850  t1 = t; /* new (best so far) W */
851  t2 = t+(1.0-t)*LAMBDA; /* new cutting plane */
852  /* t2 = t+(1.0-t)/10.0; new cutting plane */
853 
854  /* update W to be the best so far solution */
855  sq_norm_W = update_W( t1, user_data );
856 
857  /* select a new cut */
858  xi = 0;
859  cut_length = 0;
860  ocas.trn_err = 0;
861  new_b = 0;
862  for(i=0; i < nData; i++ ) {
863 
864 /* if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 ) */
865  if( (old_output[i]*(1-t2) + t2*output[i]) <= C[i] )
866  {
867  new_cut[cut_length] = i;
868  cut_length++;
869  new_b += C[i];
870  }
871 
872  output[i] = old_output[i]*(1-t1) + t1*output[i];
873 
874 /* if( output[i] <= 1) xi += 1-output[i];*/
875  if( output[i] <= C[i]) xi += C[i]-output[i];
876  if( output[i] <= 0) ocas.trn_err++;
877 
878  }
879 
880 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/
881  ocas.Q_P = 0.5*sq_norm_W + xi;
882 
883  ocas.ocas_time = get_time() - ocas_start_time;
884 
885  /* 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",
886  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),
887  ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag );
888  */
889 
890  start_time = get_time();
891  ocas_print(ocas);
892  ocas.print_time += get_time() - start_time;
893 
894  break;
895  }
896 
897  /* Stopping conditions */
898  if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
899  if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
900  if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
901  if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
902  if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
903 
904  } /* end of the main loop */
905 
906 cleanup:
907 
908  LIBOCAS_FREE(H);
909  LIBOCAS_FREE(b);
910  LIBOCAS_FREE(alpha);
911  LIBOCAS_FREE(new_cut);
912  LIBOCAS_FREE(I);
913  LIBOCAS_FREE(diag_H);
914  LIBOCAS_FREE(output);
915  LIBOCAS_FREE(old_output);
916  LIBOCAS_FREE(hpf);
917 /* LIBOCAS_FREE(hpi);*/
918  LIBOCAS_FREE(hpb);
919  LIBOCAS_FREE(Ci);
920  LIBOCAS_FREE(Bi);
921 
922  ocas.ocas_time = get_time() - ocas_start_time;
923 
924  return(ocas);
925 }
926 
927 
928 
929 /*----------------------------------------------------------------------
930  Multiclass SVM-Ocas solver
931  ----------------------------------------------------------------------*/
932 
933 /* Helper function needed by the multi-class SVM linesearch.
934 
935  - This function finds a simplified representation of a piece-wise linear function
936  by splitting the domain into intervals and fining active terms for these intevals */
937 static void findactive(float64_t *Theta, float64_t *SortedA, uint32_t *nSortedA, float64_t *A, float64_t *B, int n,
938  int (*sort)(float64_t*, float64_t*, uint32_t))
939 {
940  float64_t tmp, theta;
941  int32_t i, j, idx, idx2 = 0, start;
942 
943  sort(A,B,n);
944 
945  tmp = B[0];
946  idx = 0;
947  i = 0;
948  while( i < n-1 && A[i] == A[i+1])
949  {
950  if( B[i+1] > B[idx] )
951  {
952  idx = i+1;
953  tmp = B[i+1];
954  }
955  i++;
956  }
957 
958  (*nSortedA) = 1;
959  SortedA[0] = A[idx];
960 
961  while(1)
962  {
963  start = idx + 1;
964  while( start < n && A[idx] == A[start])
965  start++;
966 
967  theta = LIBOCAS_PLUS_INF;
968  for(j=start; j < n; j++)
969  {
970  tmp = (B[j] - B[idx])/(A[idx]-A[j]);
971  if( tmp < theta)
972  {
973  theta = tmp;
974  idx2 = j;
975  }
976  }
977 
978  if( theta < LIBOCAS_PLUS_INF)
979  {
980  Theta[(*nSortedA) - 1] = theta;
981  SortedA[(*nSortedA)] = A[idx2];
982  (*nSortedA)++;
983  idx = idx2;
984  }
985  else
986  return;
987  }
988 }
989 
990 
991 /*----------------------------------------------------------------------
992  Multiclass linear OCAS-SVM solver.
993  ----------------------------------------------------------------------*/
994 ocas_return_value_T msvm_ocas_solver(
995  float64_t C,
996  float64_t *data_y,
997  uint32_t nY,
998  uint32_t nData,
999  float64_t TolRel,
1000  float64_t TolAbs,
1001  float64_t QPBound,
1002  float64_t MaxTime,
1003  uint32_t _BufSize,
1004  uint8_t Method,
1005  void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*),
1006  float64_t (*update_W)(float64_t, void*),
1007  int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, void*),
1008  int (*compute_output)(float64_t*, void* ),
1009  int (*sort)(float64_t*, float64_t*, uint32_t),
1010  void (*ocas_print)(ocas_return_value_T),
1011  void* user_data)
1012 {
1013  ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1014  float64_t *b, *alpha, *diag_H;
1015  float64_t *output, *old_output;
1016  float64_t xi, sq_norm_W, QPSolverTolRel, QPSolverTolAbs, dot_prod_WoldW, sq_norm_oldW;
1017  float64_t A0, B0, t, t1, t2, R, tmp, element_b, x;
1018  float64_t *A, *B, *theta, *Theta, *sortedA, *Add;
1019  float64_t start_time, ocas_start_time, grad_sum, grad, min_x = 0, old_x, old_grad;
1020  uint32_t i, y, y2, ypred = 0, *new_cut, cnt1, cnt2, j, nSortedA, idx;
1021  uint32_t *I;
1022  uint8_t S = 1;
1023  libqp_state_T qp_exitflag;
1024 
1025  ocas_start_time = get_time();
1026  ocas.qp_solver_time = 0;
1027  ocas.output_time = 0;
1028  ocas.sort_time = 0;
1029  ocas.add_time = 0;
1030  ocas.w_time = 0;
1031  ocas.print_time = 0;
1032 
1033  BufSize = _BufSize;
1034 
1035  QPSolverTolRel = TolRel*0.5;
1036  QPSolverTolAbs = TolAbs*0.5;
1037 
1038  H=NULL;
1039  b=NULL;
1040  alpha=NULL;
1041  new_cut=NULL;
1042  I=NULL;
1043  diag_H=NULL;
1044  output=NULL;
1045  old_output=NULL;
1046  A = NULL;
1047  B = NULL;
1048  theta = NULL;
1049  Theta = NULL;
1050  sortedA = NULL;
1051  Add = NULL;
1052 
1053  /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
1054  H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t));
1055  if(H == NULL)
1056  {
1057  ocas.exitflag=-2;
1058  goto cleanup;
1059  }
1060 
1061  /* bias of cutting planes */
1062  b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1063  if(b == NULL)
1064  {
1065  ocas.exitflag=-2;
1066  goto cleanup;
1067  }
1068 
1069  alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1070  if(alpha == NULL)
1071  {
1072  ocas.exitflag=-2;
1073  goto cleanup;
1074  }
1075 
1076  /* indices of examples which define a new cut */
1077  new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
1078  if(new_cut == NULL)
1079  {
1080  ocas.exitflag=-2;
1081  goto cleanup;
1082  }
1083 
1084  I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
1085  if(I == NULL)
1086  {
1087  ocas.exitflag=-2;
1088  goto cleanup;
1089  }
1090 
1091  for(i=0; i< BufSize; i++)
1092  I[i] = 1;
1093 
1094  diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t));
1095  if(diag_H == NULL)
1096  {
1097  ocas.exitflag=-2;
1098  goto cleanup;
1099  }
1100 
1101  output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1102  if(output == NULL)
1103  {
1104  ocas.exitflag=-2;
1105  goto cleanup;
1106  }
1107 
1108  old_output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1109  if(old_output == NULL)
1110  {
1111  ocas.exitflag=-2;
1112  goto cleanup;
1113  }
1114 
1115  /* auxciliary variables used in the linesearch */
1116  A = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1117  if(A == NULL)
1118  {
1119  ocas.exitflag=-2;
1120  goto cleanup;
1121  }
1122 
1123  B = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1124  if(B == NULL)
1125  {
1126  ocas.exitflag=-2;
1127  goto cleanup;
1128  }
1129 
1130  theta = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t));
1131  if(theta == NULL)
1132  {
1133  ocas.exitflag=-2;
1134  goto cleanup;
1135  }
1136 
1137  sortedA = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t));
1138  if(sortedA == NULL)
1139  {
1140  ocas.exitflag=-2;
1141  goto cleanup;
1142  }
1143 
1144  Theta = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1145  if(Theta == NULL)
1146  {
1147  ocas.exitflag=-2;
1148  goto cleanup;
1149  }
1150 
1151  Add = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t));
1152  if(Add == NULL)
1153  {
1154  ocas.exitflag=-2;
1155  goto cleanup;
1156  }
1157 
1158  /* Set initial values*/
1159  ocas.nCutPlanes = 0;
1160  ocas.exitflag = 0;
1161  ocas.nIter = 0;
1162  ocas.Q_D = 0;
1163  ocas.trn_err = nData;
1164  R = (float64_t)nData;
1165  sq_norm_W = 0;
1166  element_b = (float64_t)nData;
1167  ocas.Q_P = 0.5*sq_norm_W + C*R;
1168 
1169  /* initial cutting plane */
1170  for(i=0; i < nData; i++)
1171  {
1172  y2 = (uint32_t)data_y[i]-1;
1173 
1174  if(y2 > 0)
1175  new_cut[i] = 0;
1176  else
1177  new_cut[i] = 1;
1178 
1179  }
1180 
1181  ocas.ocas_time = get_time() - ocas_start_time;
1182 
1183  start_time = get_time();
1184  ocas_print(ocas);
1185  ocas.print_time += get_time() - start_time;
1186 
1187  /* main loop of the OCAS */
1188  while( ocas.exitflag == 0 )
1189  {
1190  ocas.nIter++;
1191 
1192  /* append a new cut to the buffer and update H */
1193  b[ocas.nCutPlanes] = -(float64_t)element_b;
1194 
1195  start_time = get_time();
1196 
1197  if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, ocas.nCutPlanes, user_data ) != 0)
1198  {
1199  ocas.exitflag=-2;
1200  goto cleanup;
1201  }
1202 
1203  ocas.add_time += get_time() - start_time;
1204 
1205  /* copy newly appended row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
1206  diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
1207  for(i=0; i < ocas.nCutPlanes; i++)
1208  {
1209  H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
1210  }
1211 
1212  ocas.nCutPlanes++;
1213 
1214  /* call inner QP solver */
1215  start_time = get_time();
1216 
1217  qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
1218  ocas.nCutPlanes, QPSolverMaxIter, QPSolverTolAbs, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
1219 
1220  ocas.qp_exitflag = qp_exitflag.exitflag;
1221 
1222  ocas.qp_solver_time += get_time() - start_time;
1223  ocas.Q_D = -qp_exitflag.QP;
1224 
1225  ocas.nNZAlpha = 0;
1226  for(i=0; i < ocas.nCutPlanes; i++)
1227  if( alpha[i] != 0) ocas.nNZAlpha++;
1228 
1229  sq_norm_oldW = sq_norm_W;
1230  start_time = get_time();
1231  compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
1232  ocas.w_time += get_time() - start_time;
1233 
1234  /* select a new cut */
1235  switch( Method )
1236  {
1237  /* cutting plane algorithm implemented in SVMperf and BMRM */
1238  case 0:
1239 
1240  start_time = get_time();
1241  if( compute_output( output, user_data ) != 0)
1242  {
1243  ocas.exitflag=-2;
1244  goto cleanup;
1245  }
1246  ocas.output_time += get_time()-start_time;
1247 
1248  /* the following loop computes: */
1249  element_b = 0.0; /* element_b = R(old_W) - g'*old_W */
1250  R = 0; /* R(W) = sum_i max_y ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */
1251  ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */
1252  /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */
1253  for(i=0; i < nData; i++)
1254  {
1255  y2 = (uint32_t)data_y[i]-1;
1256 
1257  for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1258  {
1259  if(y2 != y && xi < output[LIBOCAS_INDEX(y,i,nY)])
1260  {
1261  xi = output[LIBOCAS_INDEX(y,i,nY)];
1262  ypred = y;
1263  }
1264  }
1265 
1266  if(xi >= output[LIBOCAS_INDEX(y2,i,nY)])
1267  ocas.trn_err ++;
1268 
1269  xi = LIBOCAS_MAX(0,xi+1-output[LIBOCAS_INDEX(y2,i,nY)]);
1270  R += xi;
1271  if(xi > 0)
1272  {
1273  element_b++;
1274  new_cut[i] = ypred;
1275  }
1276  else
1277  new_cut[i] = y2;
1278  }
1279 
1280  ocas.Q_P = 0.5*sq_norm_W + C*R;
1281 
1282  ocas.ocas_time = get_time() - ocas_start_time;
1283 
1284  start_time = get_time();
1285  ocas_print(ocas);
1286  ocas.print_time += get_time() - start_time;
1287 
1288  break;
1289 
1290  /* The OCAS solver */
1291  case 1:
1292  memcpy( old_output, output, sizeof(float64_t)*nData*nY );
1293 
1294  start_time = get_time();
1295  if( compute_output( output, user_data ) != 0)
1296  {
1297  ocas.exitflag=-2;
1298  goto cleanup;
1299  }
1300  ocas.output_time += get_time()-start_time;
1301 
1302  A0 = sq_norm_W - 2*dot_prod_WoldW + sq_norm_oldW;
1303  B0 = dot_prod_WoldW - sq_norm_oldW;
1304 
1305  for(i=0; i < nData; i++)
1306  {
1307  y2 = (uint32_t)data_y[i]-1;
1308 
1309  for(y=0; y < nY; y++)
1310  {
1311  A[LIBOCAS_INDEX(y,i,nY)] = C*(output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y,i,nY)]
1312  + old_output[LIBOCAS_INDEX(y2,i,nY)] - output[LIBOCAS_INDEX(y2,i,nY)]);
1313  B[LIBOCAS_INDEX(y,i,nY)] = C*(old_output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y2,i,nY)]
1314  + (float64_t)(y != y2));
1315  }
1316  }
1317 
1318  /* linesearch */
1319 /* new_x = msvm_linesearch_mex(A0,B0,AA*C,BB*C);*/
1320 
1321  grad_sum = B0;
1322  cnt1 = 0;
1323  cnt2 = 0;
1324  for(i=0; i < nData; i++)
1325  {
1326  findactive(theta,sortedA,&nSortedA,&A[i*nY],&B[i*nY],nY,sort);
1327 
1328  idx = 0;
1329  while( idx < nSortedA-1 && theta[idx] < 0 )
1330  idx++;
1331 
1332  grad_sum += sortedA[idx];
1333 
1334  for(j=idx; j < nSortedA-1; j++)
1335  {
1336  Theta[cnt1] = theta[j];
1337  cnt1++;
1338  }
1339 
1340  for(j=idx+1; j < nSortedA; j++)
1341  {
1342  Add[cnt2] = -sortedA[j-1]+sortedA[j];
1343  cnt2++;
1344  }
1345  }
1346 
1347  start_time = get_time();
1348  sort(Theta,Add,cnt1);
1349  ocas.sort_time += get_time() - start_time;
1350 
1351  grad = grad_sum;
1352  if(grad >= 0)
1353  {
1354  min_x = 0;
1355  }
1356  else
1357  {
1358  old_x = 0;
1359  old_grad = grad;
1360 
1361  for(i=0; i < cnt1; i++)
1362  {
1363  x = Theta[i];
1364 
1365  grad = x*A0 + grad_sum;
1366 
1367  if(grad >=0)
1368  {
1369 
1370  min_x = (grad*old_x - old_grad*x)/(grad - old_grad);
1371 
1372  break;
1373  }
1374  else
1375  {
1376  grad_sum = grad_sum + Add[i];
1377 
1378  grad = x*A0 + grad_sum;
1379  if( grad >= 0)
1380  {
1381  min_x = x;
1382  break;
1383  }
1384  }
1385 
1386  old_grad = grad;
1387  old_x = x;
1388  }
1389  }
1390  /* end of the linesearch which outputs min_x */
1391 
1392  t = min_x;
1393  t1 = t; /* new (best so far) W */
1394  t2 = t+(1.0-t)*LAMBDA; /* new cutting plane */
1395  /* t2 = t+(1.0-t)/10.0; */
1396 
1397  /* update W to be the best so far solution */
1398  sq_norm_W = update_W( t1, user_data );
1399 
1400  /* the following code computes a new cutting plane: */
1401  element_b = 0.0; /* element_b = R(old_W) - g'*old_W */
1402  /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */
1403  for(i=0; i < nData; i++)
1404  {
1405  y2 = (uint32_t)data_y[i]-1;
1406 
1407  for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1408  {
1409  tmp = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y,i,nY)];
1410  if(y2 != y && xi < tmp)
1411  {
1412  xi = tmp;
1413  ypred = y;
1414  }
1415  }
1416 
1417  tmp = old_output[LIBOCAS_INDEX(y2,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y2,i,nY)];
1418  xi = LIBOCAS_MAX(0,xi+1-tmp);
1419  if(xi > 0)
1420  {
1421  element_b++;
1422  new_cut[i] = ypred;
1423  }
1424  else
1425  new_cut[i] = y2;
1426  }
1427 
1428  /* compute Risk, class. error and update outputs to correspond to the new W */
1429  ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */
1430  R = 0;
1431  for(i=0; i < nData; i++)
1432  {
1433  y2 = (uint32_t)data_y[i]-1;
1434 
1435  for(tmp=-LIBOCAS_PLUS_INF, y=0; y < nY; y++)
1436  {
1437  output[LIBOCAS_INDEX(y,i,nY)] = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t1) + t1*output[LIBOCAS_INDEX(y,i,nY)];
1438 
1439  if(y2 != y && tmp < output[LIBOCAS_INDEX(y,i,nY)])
1440  {
1441  ypred = y;
1442  tmp = output[LIBOCAS_INDEX(y,i,nY)];
1443  }
1444  }
1445 
1446  R += LIBOCAS_MAX(0,1+tmp - output[LIBOCAS_INDEX(y2,i,nY)]);
1447  if( tmp >= output[LIBOCAS_INDEX(y2,i,nY)])
1448  ocas.trn_err ++;
1449  }
1450 
1451  ocas.Q_P = 0.5*sq_norm_W + C*R;
1452 
1453 
1454  /* get time and print status */
1455  ocas.ocas_time = get_time() - ocas_start_time;
1456 
1457  start_time = get_time();
1458  ocas_print(ocas);
1459  ocas.print_time += get_time() - start_time;
1460 
1461  break;
1462 
1463  }
1464 
1465  /* Stopping conditions */
1466  if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
1467  if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
1468  if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
1469  if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
1470  if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
1471 
1472  } /* end of the main loop */
1473 
1474 cleanup:
1475 
1476  LIBOCAS_FREE(H);
1477  LIBOCAS_FREE(b);
1478  LIBOCAS_FREE(alpha);
1479  LIBOCAS_FREE(new_cut);
1480  LIBOCAS_FREE(I);
1481  LIBOCAS_FREE(diag_H);
1482  LIBOCAS_FREE(output);
1483  LIBOCAS_FREE(old_output);
1484  LIBOCAS_FREE(A);
1485  LIBOCAS_FREE(B);
1486  LIBOCAS_FREE(theta);
1487  LIBOCAS_FREE(Theta);
1488  LIBOCAS_FREE(sortedA);
1489  LIBOCAS_FREE(Add);
1490 
1491  ocas.ocas_time = get_time() - ocas_start_time;
1492 
1493  return(ocas);
1494 }
1495 }

SHOGUN Machine Learning Toolbox - Documentation