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