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