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