SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
libppbm.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  * libp3bm.h: Implementation of the Proximal Point P-BMRM solver for SO training
8  *
9  * Copyright (C) 2012 Michal Uricar, uricamic@cmp.felk.cvut.cz
10  *
11  * Implementation of the Proximal Point P-BMRM
12  *--------------------------------------------------------------------- */
13 
15 #include <shogun/lib/external/libqp.h>
16 #include <shogun/lib/Time.h>
17 
18 namespace shogun
19 {
20 static const uint32_t QPSolverMaxIter=0xFFFFFFFF;
21 static const float64_t epsilon=0.0;
22 
23 static float64_t *H, *H2;
24 
25 /*----------------------------------------------------------------------
26  Returns pointer at i-th column of Hessian matrix.
27  ----------------------------------------------------------------------*/
28 static const float64_t *get_col( uint32_t i)
29 {
30  return( &H2[ BufSize*i ] );
31 }
32 
34  CDualLibQPBMSOSVM *machine,
35  float64_t* W,
36  float64_t TolRel,
37  float64_t TolAbs,
38  float64_t _lambda,
39  uint32_t _BufSize,
40  bool cleanICP,
41  uint32_t cleanAfter,
42  float64_t K,
43  uint32_t Tmax,
44  bool verbose)
45 {
46  BmrmStatistics ppbmrm;
47  libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0};
48  float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2;
49  float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0;
50  float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0;
51  float64_t lastFp, wdist, gamma=0.0;
52  floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps;
53  uint32_t *I, *I2, *I_start, *I_good;
54  uint8_t S=1;
55  CStructuredModel* model=machine->get_model();
56  uint32_t nDim=model->get_dim();
57  CSOSVMHelper* helper = NULL;
58  uint32_t qp_cnt=0;
59  bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
60  float64_t *A_1=NULL, *A_2=NULL;
61  bool *map=NULL, tuneAlpha=true, flag=true, alphaChanged=false, isThereGoodSolution=false;
62 
63  CTime ttime;
64  float64_t tstart, tstop;
65 
66 
67  tstart=ttime.cur_time_diff(false);
68 
69  BufSize=_BufSize;
70  QPSolverTolRel=1e-9;
71 
72  H=NULL;
73  b=NULL;
74  beta=NULL;
75  A=NULL;
76  subgrad=NULL;
77  diag_H=NULL;
78  I=NULL;
79  prevW=NULL;
80  wt=NULL;
81  diag_H2=NULL;
82  b2=NULL;
83  I2=NULL;
84  H2=NULL;
85  I_good=NULL;
86  I_start=NULL;
87  beta_start=NULL;
88  beta_good=NULL;
89 
90  alpha=0.0;
91 
93 
94  ASSERT(nDim > 0);
95  ASSERT(BufSize > 0);
96  REQUIRE(BufSize < (std::numeric_limits<size_t>::max() / nDim),
97  "overflow: %u * %u > %u -- biggest possible BufSize=%u or nDim=%u\n",
98  BufSize, nDim, std::numeric_limits<size_t>::max(),
99  (std::numeric_limits<size_t>::max() / nDim),
100  (std::numeric_limits<size_t>::max() / BufSize));
101 
102  A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, float64_t);
103 
104  b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
105 
106  beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
107 
108  subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
109 
110  diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
111 
112  I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
113 
114  /* structure for maintaining inactive CPs info */
115  ICP_stats icp_stats;
116  icp_stats.maxCPs = BufSize;
117  icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
118  icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*);
119  icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
120 
121  cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
122 
123  prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
124 
125  wt= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
126 
127  if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL ||
128  diag_H==NULL || I==NULL || icp_stats.ICPcounter==NULL ||
129  icp_stats.ICPs==NULL || icp_stats.ACPs==NULL ||
130  cp_list==NULL || prevW==NULL || wt==NULL)
131  {
132  ppbmrm.exitflag=-2;
133  goto cleanup;
134  }
135 
136  map= (bool*) LIBBMRM_CALLOC(BufSize, bool);
137 
138  if (map==NULL)
139  {
140  ppbmrm.exitflag=-2;
141  goto cleanup;
142  }
143 
144  memset( (bool*) map, true, BufSize);
145 
146  /* Temporary buffers for ICP removal */
147  icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
148 
149  if (icp_stats.H_buff==NULL)
150  {
151  ppbmrm.exitflag=-2;
152  goto cleanup;
153  }
154 
155  /* Temporary buffers */
156  beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
157 
158  beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
159 
160  b2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
161 
162  diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
163 
164  H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
165 
166  I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
167 
168  I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
169 
170  I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
171 
172  if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
173  I_start==NULL || I_good==NULL || I2==NULL || H2==NULL)
174  {
175  ppbmrm.exitflag=-2;
176  goto cleanup;
177  }
178 
179  ppbmrm.hist_Fp.resize_vector(BufSize);
180  ppbmrm.hist_Fd.resize_vector(BufSize);
181  ppbmrm.hist_wdist.resize_vector(BufSize);
182 
183  /* Iinitial solution */
184  R = machine->risk(subgrad, W);
185 
186  ppbmrm.nCP=0;
187  ppbmrm.nIter=0;
188  ppbmrm.exitflag=0;
189 
190  b[0]=-R;
191 
192  /* Cutting plane auxiliary double linked list */
193  LIBBMRM_MEMCPY(A, subgrad, nDim*sizeof(float64_t));
194  map[0]=false;
195  cp_list->address=&A[0];
196  cp_list->idx=0;
197  cp_list->prev=NULL;
198  cp_list->next=NULL;
199  CPList_head=cp_list;
200  CPList_tail=cp_list;
201 
202  /* Compute initial value of Fp, Fd, assuming that W is zero vector */
203  sq_norm_Wdiff=0.0;
204 
205  b[0] = SGVector<float64_t>::dot(subgrad, W, nDim);
206  sq_norm_W = SGVector<float64_t>::dot(W, W, nDim);
207  for (uint32_t j=0; j<nDim; ++j)
208  {
209  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
210  }
211 
212  ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
213  ppbmrm.Fd=-LIBBMRM_PLUS_INF;
214  lastFp=ppbmrm.Fp;
215  wdist=CMath::sqrt(sq_norm_Wdiff);
216 
217  K = (sq_norm_W == 0.0) ? 0.4 : 0.01*CMath::sqrt(sq_norm_W);
218 
219  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
220 
221  tstop=ttime.cur_time_diff(false);
222 
223  /* Keep history of Fp, Fd, wdist */
224  ppbmrm.hist_Fp[0]=ppbmrm.Fp;
225  ppbmrm.hist_Fd[0]=ppbmrm.Fd;
226  ppbmrm.hist_wdist[0]=wdist;
227 
228  /* Verbose output */
229 
230  if (verbose)
231  SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n",
232  ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, R, K);
233 
234  if (verbose)
235  helper = machine->get_helper();
236 
237  /* main loop */
238 
239  while (ppbmrm.exitflag==0)
240  {
241  tstart=ttime.cur_time_diff(false);
242  ppbmrm.nIter++;
243 
244  /* Update H */
245 
246  if (ppbmrm.nCP>0)
247  {
248  A_2=get_cutting_plane(CPList_tail);
249  cp_ptr=CPList_head;
250 
251  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
252  {
253  A_1=get_cutting_plane(cp_ptr);
254  cp_ptr=cp_ptr->next;
255  rsum=SGVector<float64_t>::dot(A_1, A_2, nDim);
256 
257  H[LIBBMRM_INDEX(ppbmrm.nCP, i, BufSize)]
258  = H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)]
259  = rsum;
260  }
261  }
262 
263  A_2=get_cutting_plane(CPList_tail);
264  rsum = SGVector<float64_t>::dot(A_2, A_2, nDim);
265 
266  H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]=rsum;
267 
268  diag_H[ppbmrm.nCP]=H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)];
269  I[ppbmrm.nCP]=1;
270 
271  beta[ppbmrm.nCP]=0.0; // [beta; 0]
272  ppbmrm.nCP++;
273 
274  /* tune alpha cycle */
275  /* ---------------------------------------------------------------------- */
276 
277  flag=true;
278  isThereGoodSolution=false;
279  LIBBMRM_MEMCPY(beta_start, beta, ppbmrm.nCP*sizeof(float64_t));
280  LIBBMRM_MEMCPY(I_start, I, ppbmrm.nCP*sizeof(uint32_t));
281  qp_cnt=0;
282  alpha_good=alpha;
283 
284  if (tuneAlpha)
285  {
286  alpha_start=alpha; alpha=0.0;
287  beta[ppbmrm.nCP]=0.0;
288  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
289  I2[ppbmrm.nCP]=1;
290 
291  /* add alpha-dependent terms to H, diag_h and b */
292  cp_ptr=CPList_head;
293 
294  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
295  {
296  A_1=get_cutting_plane(cp_ptr);
297  cp_ptr=cp_ptr->next;
298 
299  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
300 
301  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
302  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
303 
304  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
305  H2[LIBBMRM_INDEX(i, j, BufSize)]=
306  H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
307  }
308 
309  /* solve QP with current alpha */
310  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
311  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
312  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
313  qp_cnt++;
314  Fd_alpha0=-qp_exitflag.QP;
315 
316  /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */
317  for (uint32_t i=0; i<nDim; ++i)
318  {
319  rsum=0.0;
320  cp_ptr=CPList_head;
321 
322  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
323  {
324  A_1=get_cutting_plane(cp_ptr);
325  cp_ptr=cp_ptr->next;
326  rsum+=A_1[i]*beta[j];
327  }
328 
329  wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
330  }
331 
332  sq_norm_Wdiff=0.0;
333 
334  for (uint32_t i=0; i<nDim; ++i)
335  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
336 
337  if (CMath::sqrt(sq_norm_Wdiff) <= K)
338  {
339  flag=false;
340 
341  if (alpha!=alpha_start)
342  alphaChanged=true;
343  }
344  else
345  {
346  alpha=alpha_start;
347  }
348 
349  while(flag)
350  {
351  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
352  LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
353  I2[ppbmrm.nCP]=1;
354  beta[ppbmrm.nCP]=0.0;
355 
356  /* add alpha-dependent terms to H, diag_h and b */
357  cp_ptr=CPList_head;
358 
359  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
360  {
361  A_1=get_cutting_plane(cp_ptr);
362  cp_ptr=cp_ptr->next;
363 
364  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
365 
366  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
367  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
368 
369  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
370  H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
371  }
372 
373  /* solve QP with current alpha */
374  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
375  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
376  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
377  qp_cnt++;
378 
379  /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */
380  for (uint32_t i=0; i<nDim; ++i)
381  {
382  rsum=0.0;
383  cp_ptr=CPList_head;
384 
385  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
386  {
387  A_1=get_cutting_plane(cp_ptr);
388  cp_ptr=cp_ptr->next;
389  rsum+=A_1[i]*beta[j];
390  }
391 
392  wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
393  }
394 
395  sq_norm_Wdiff=0.0;
396  for (uint32_t i=0; i<nDim; ++i)
397  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
398 
399  if (CMath::sqrt(sq_norm_Wdiff) > K)
400  {
401  /* if there is a record of some good solution
402  * (i.e. adjust alpha by division by 2) */
403 
404  if (isThereGoodSolution)
405  {
406  LIBBMRM_MEMCPY(beta, beta_good, ppbmrm.nCP*sizeof(float64_t));
407  LIBBMRM_MEMCPY(I2, I_good, ppbmrm.nCP*sizeof(uint32_t));
408  alpha=alpha_good;
409  qp_exitflag=qp_exitflag_good;
410  flag=false;
411  }
412  else
413  {
414  if (alpha == 0)
415  {
416  alpha=1.0;
417  alphaChanged=true;
418  }
419  else
420  {
421  alpha*=2;
422  alphaChanged=true;
423  }
424  }
425  }
426  else
427  {
428  if (alpha > 0)
429  {
430  /* keep good solution and try for alpha /= 2 if previous alpha was 1 */
431  LIBBMRM_MEMCPY(beta_good, beta, ppbmrm.nCP*sizeof(float64_t));
432  LIBBMRM_MEMCPY(I_good, I2, ppbmrm.nCP*sizeof(uint32_t));
433  alpha_good=alpha;
434  qp_exitflag_good=qp_exitflag;
435  isThereGoodSolution=true;
436 
437  if (alpha!=1.0)
438  {
439  alpha/=2.0;
440  alphaChanged=true;
441  }
442  else
443  {
444  alpha=0.0;
445  alphaChanged=true;
446  }
447  }
448  else
449  {
450  flag=false;
451  }
452  }
453  }
454  }
455  else
456  {
457  alphaChanged=false;
458  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
459  LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
460 
461  /* add alpha-dependent terms to H, diag_h and b */
462  cp_ptr=CPList_head;
463 
464  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
465  {
466  A_1=get_cutting_plane(cp_ptr);
467  cp_ptr=cp_ptr->next;
468 
469  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
470 
471  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
472  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
473 
474  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
475  H2[LIBBMRM_INDEX(i, j, BufSize)]=
476  H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
477  }
478  /* solve QP with current alpha */
479  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
480  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
481  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
482  qp_cnt++;
483  }
484 
485  /* ----------------------------------------------------------------------------------------------- */
486 
487  /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */
488  ppbmrm.nzA=0;
489 
490  for (uint32_t aaa=0; aaa<ppbmrm.nCP; ++aaa)
491  {
492  if (beta[aaa]>epsilon)
493  {
494  ++ppbmrm.nzA;
495  icp_stats.ICPcounter[aaa]=0;
496  }
497  else
498  {
499  icp_stats.ICPcounter[aaa]+=1;
500  }
501  }
502 
503  /* W update */
504  for (uint32_t i=0; i<nDim; ++i)
505  {
506  rsum=0.0;
507  cp_ptr=CPList_head;
508 
509  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
510  {
511  A_1=get_cutting_plane(cp_ptr);
512  cp_ptr=cp_ptr->next;
513  rsum+=A_1[i]*beta[j];
514  }
515 
516  W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha);
517  }
518 
519  /* risk and subgradient computation */
520  R = machine->risk(subgrad, W);
521  add_cutting_plane(&CPList_tail, map, A,
522  find_free_idx(map, BufSize), subgrad, nDim);
523 
524  sq_norm_W=SGVector<float64_t>::dot(W, W, nDim);
525  sq_norm_prevW=SGVector<float64_t>::dot(prevW, prevW, nDim);
526  b[ppbmrm.nCP]=SGVector<float64_t>::dot(subgrad, W, nDim) - R;
527 
528  sq_norm_Wdiff=0.0;
529  for (uint32_t j=0; j<nDim; ++j)
530  {
531  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
532  }
533 
534  /* compute Fp and Fd */
535  ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
536  ppbmrm.Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
537 
538  /* gamma + tuneAlpha flag */
539  if (alphaChanged)
540  {
541  eps=1.0-(ppbmrm.Fd/ppbmrm.Fp);
542  gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
543  }
544 
545  if ((lastFp-ppbmrm.Fp) <= gamma)
546  {
547  tuneAlpha=true;
548  }
549  else
550  {
551  tuneAlpha=false;
552  }
553 
554  /* Stopping conditions - set only with nonzero alpha */
555  if (alpha==0.0)
556  {
557  if (ppbmrm.Fp-ppbmrm.Fd<=TolRel*LIBBMRM_ABS(ppbmrm.Fp))
558  ppbmrm.exitflag=1;
559 
560  if (ppbmrm.Fp-ppbmrm.Fd<=TolAbs)
561  ppbmrm.exitflag=2;
562  }
563 
564  if (ppbmrm.nCP>=BufSize)
565  ppbmrm.exitflag=-1;
566 
567  tstop=ttime.cur_time_diff(false);
568 
569  /* compute wdist (= || W_{t+1} - W_{t} || ) */
570  sq_norm_Wdiff=0.0;
571 
572  for (uint32_t i=0; i<nDim; ++i)
573  {
574  sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
575  }
576 
577  wdist=CMath::sqrt(sq_norm_Wdiff);
578 
579  /* Keep history of Fp, Fd, wdist */
580  ppbmrm.hist_Fp[ppbmrm.nIter]=ppbmrm.Fp;
581  ppbmrm.hist_Fd[ppbmrm.nIter]=ppbmrm.Fd;
582  ppbmrm.hist_wdist[ppbmrm.nIter]=wdist;
583 
584  /* Verbose output */
585  if (verbose)
586  SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, (Fp-Fd)=%lf, (Fp-Fd)/Fp=%lf, R=%lf, nCP=%d, nzA=%d, wdist=%lf, alpha=%lf, qp_cnt=%d, gamma=%lf, tuneAlpha=%d\n",
587  ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, ppbmrm.Fp-ppbmrm.Fd,
588  (ppbmrm.Fp-ppbmrm.Fd)/ppbmrm.Fp, R, ppbmrm.nCP, ppbmrm.nzA, wdist, alpha,
589  qp_cnt, gamma, tuneAlpha);
590 
591  /* Check size of Buffer */
592  if (ppbmrm.nCP>=BufSize)
593  {
594  ppbmrm.exitflag=-2;
595  SG_SERROR("Buffer exceeded.\n")
596  }
597 
598  /* keep w_t + Fp */
599  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
600  lastFp=ppbmrm.Fp;
601 
602  /* Inactive Cutting Planes (ICP) removal */
603  if (cleanICP)
604  {
605  clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
606  }
607 
608  // next CP would exceed BufSize
609  if (ppbmrm.nCP+1 >= BufSize)
610  ppbmrm.exitflag=-1;
611 
612  /* Debug: compute objective and training error */
613  if (verbose)
614  {
615  SGVector<float64_t> w_debug(W, nDim, false);
616  float64_t primal = CSOSVMHelper::primal_objective(w_debug, model, _lambda);
617  float64_t train_error = CSOSVMHelper::average_loss(w_debug, model);
618  helper->add_debug_info(primal, ppbmrm.nIter, train_error);
619  }
620  } /* end of main loop */
621 
622  if (verbose)
623  {
624  helper->terminate();
625  SG_UNREF(helper);
626  }
627 
628  ppbmrm.hist_Fp.resize_vector(ppbmrm.nIter);
629  ppbmrm.hist_Fd.resize_vector(ppbmrm.nIter);
630  ppbmrm.hist_wdist.resize_vector(ppbmrm.nIter);
631 
632  cp_ptr=CPList_head;
633 
634  while(cp_ptr!=NULL)
635  {
636  cp_ptr2=cp_ptr;
637  cp_ptr=cp_ptr->next;
638  LIBBMRM_FREE(cp_ptr2);
639  cp_ptr2=NULL;
640  }
641 
642  cp_list=NULL;
643 
644 cleanup:
645 
646  LIBBMRM_FREE(H);
647  LIBBMRM_FREE(b);
648  LIBBMRM_FREE(beta);
649  LIBBMRM_FREE(A);
650  LIBBMRM_FREE(subgrad);
651  LIBBMRM_FREE(diag_H);
652  LIBBMRM_FREE(I);
653  LIBBMRM_FREE(icp_stats.ICPcounter);
654  LIBBMRM_FREE(icp_stats.ICPs);
655  LIBBMRM_FREE(icp_stats.ACPs);
656  LIBBMRM_FREE(icp_stats.H_buff);
657  LIBBMRM_FREE(map);
658  LIBBMRM_FREE(prevW);
659  LIBBMRM_FREE(wt);
660  LIBBMRM_FREE(beta_start);
661  LIBBMRM_FREE(beta_good);
662  LIBBMRM_FREE(I_start);
663  LIBBMRM_FREE(I_good);
664  LIBBMRM_FREE(I2);
665  LIBBMRM_FREE(b2);
666  LIBBMRM_FREE(diag_H2);
667  LIBBMRM_FREE(H2);
668 
669  if (cp_list)
670  LIBBMRM_FREE(cp_list);
671 
672  SG_UNREF(model);
673 
674  return(ppbmrm);
675 }
676 }

SHOGUN Machine Learning Toolbox - Documentation