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

SHOGUN Machine Learning Toolbox - Documentation