SHOGUN
v2.0.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * libp3bm.h: Implementation of the Proximal Point P-BMRM solver for SO training 00008 * 00009 * Copyright (C) 2012 Michal Uricar, uricamic@cmp.felk.cvut.cz 00010 * 00011 * Implementation of the Proximal Point P-BMRM 00012 *--------------------------------------------------------------------- */ 00013 00014 #include <stdlib.h> 00015 #include <string.h> 00016 #include <math.h> 00017 00018 #include <shogun/structure/libp3bm.h> 00019 #include <shogun/lib/external/libqp.h> 00020 #include <shogun/lib/Time.h> 00021 00022 namespace shogun 00023 { 00024 static const uint32_t QPSolverMaxIter=0xFFFFFFFF; 00025 static const float64_t epsilon=0.0; 00026 00027 static float64_t *H, *H2; 00028 static uint32_t BufSize; 00029 00030 /*---------------------------------------------------------------------- 00031 Returns pointer at i-th column of Hessian matrix. 00032 ----------------------------------------------------------------------*/ 00033 static const float64_t *get_col( uint32_t i) 00034 { 00035 return( &H2[ BufSize*i ] ); 00036 } 00037 00038 bmrm_return_value_T svm_ppbm_solver( 00039 CStructuredModel* model, 00040 float64_t* W, 00041 float64_t TolRel, 00042 float64_t TolAbs, 00043 float64_t _lambda, 00044 uint32_t _BufSize, 00045 bool cleanICP, 00046 uint32_t cleanAfter, 00047 float64_t K, 00048 uint32_t Tmax, 00049 bool verbose) 00050 { 00051 bmrm_return_value_T ppbmrm; 00052 libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0}; 00053 float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2; 00054 float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0; 00055 float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0; 00056 float64_t lastFp, wdist, gamma=0.0; 00057 floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps; 00058 uint32_t *I, *I2, *I_start, *I_good, *ICPcounter, *ACPs, cntICP=0, cntACP=0; 00059 uint8_t S=1; 00060 uint32_t nDim=model->get_dim(); 00061 float64_t **ICPs; 00062 uint32_t nCP_new=0, qp_cnt=0; 00063 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL; 00064 float64_t *A_1=NULL, *A_2=NULL, *H_buff; 00065 bool *map=NULL, tuneAlpha=true, flag=true, alphaChanged=false, isThereGoodSolution=false; 00066 00067 CTime ttime; 00068 float64_t tstart, tstop; 00069 00070 00071 tstart=ttime.cur_time_diff(false); 00072 00073 BufSize=_BufSize; 00074 QPSolverTolRel=1e-9; 00075 00076 H=NULL; 00077 b=NULL; 00078 beta=NULL; 00079 A=NULL; 00080 subgrad=NULL; 00081 diag_H=NULL; 00082 I=NULL; 00083 ICPcounter=NULL; 00084 ICPs=NULL; 00085 ACPs=NULL; 00086 prevW=NULL; 00087 wt=NULL; 00088 H_buff=NULL; 00089 diag_H2=NULL; 00090 b2=NULL; 00091 I2=NULL; 00092 H2=NULL; 00093 I_good=NULL; 00094 I_start=NULL; 00095 beta_start=NULL; 00096 beta_good=NULL; 00097 00098 alpha=0.0; 00099 00100 H= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, sizeof(float64_t)); 00101 00102 A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, sizeof(float64_t)); 00103 00104 b= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00105 00106 beta= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00107 00108 subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, sizeof(float64_t)); 00109 00110 diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00111 00112 I= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00113 00114 ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00115 00116 ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, sizeof(float64_t*)); 00117 00118 ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00119 00120 cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, sizeof(bmrm_ll)); 00121 00122 prevW= (float64_t*) LIBBMRM_CALLOC(nDim, sizeof(float64_t)); 00123 00124 wt= (float64_t*) LIBBMRM_CALLOC(nDim, sizeof(float64_t)); 00125 00126 if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL || 00127 diag_H==NULL || I==NULL || ICPcounter==NULL || ICPs==NULL || 00128 ACPs==NULL || cp_list==NULL || prevW==NULL || wt==NULL) 00129 { 00130 ppbmrm.exitflag=-2; 00131 goto cleanup; 00132 } 00133 00134 map= (bool*) LIBBMRM_CALLOC(BufSize, sizeof(bool)); 00135 00136 if (map==NULL) 00137 { 00138 ppbmrm.exitflag=-2; 00139 goto cleanup; 00140 } 00141 00142 memset( (bool*) map, true, BufSize); 00143 00144 /* Temporary buffers for ICP removal */ 00145 H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, sizeof(float64_t)); 00146 00147 if (H_buff==NULL) 00148 { 00149 ppbmrm.exitflag=-2; 00150 goto cleanup; 00151 } 00152 00153 /* Temporary buffers */ 00154 beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00155 00156 beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00157 00158 b2= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00159 00160 diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00161 00162 H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, sizeof(float64_t)); 00163 00164 I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00165 00166 I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00167 00168 I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00169 00170 if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL || 00171 I_start==NULL || I_good==NULL || I2==NULL || H2==NULL) 00172 { 00173 ppbmrm.exitflag=-2; 00174 goto cleanup; 00175 } 00176 00177 ppbmrm.hist_Fp.resize_vector(BufSize); 00178 ppbmrm.hist_Fd.resize_vector(BufSize); 00179 ppbmrm.hist_wdist.resize_vector(BufSize); 00180 00181 /* Iinitial solution */ 00182 R = model->risk(subgrad, W); 00183 00184 ppbmrm.nCP=0; 00185 ppbmrm.nIter=0; 00186 ppbmrm.exitflag=0; 00187 00188 b[0]=-R; 00189 00190 /* Cutting plane auxiliary double linked list */ 00191 LIBBMRM_MEMCPY(A, subgrad, nDim*sizeof(float64_t)); 00192 map[0]=false; 00193 cp_list->address=&A[0]; 00194 cp_list->idx=0; 00195 cp_list->prev=NULL; 00196 cp_list->next=NULL; 00197 CPList_head=cp_list; 00198 CPList_tail=cp_list; 00199 00200 /* Compute initial value of Fp, Fd, assuming that W is zero vector */ 00201 sq_norm_W=0.0; 00202 sq_norm_Wdiff=0.0; 00203 00204 for (uint32_t j=0; j<nDim; ++j) 00205 { 00206 b[0]+=subgrad[j]*W[j]; 00207 sq_norm_W+=W[j]*W[j]; 00208 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]); 00209 } 00210 00211 ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff; 00212 ppbmrm.Fd=-LIBBMRM_PLUS_INF; 00213 lastFp=ppbmrm.Fp; 00214 wdist=::sqrt(sq_norm_Wdiff); 00215 00216 K = (sq_norm_W == 0.0) ? 0.4 : 0.01*::sqrt(sq_norm_W); 00217 00218 LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t)); 00219 00220 tstop=ttime.cur_time_diff(false); 00221 00222 /* Keep history of Fp, Fd, wdist */ 00223 ppbmrm.hist_Fp[0]=ppbmrm.Fp; 00224 ppbmrm.hist_Fd[0]=ppbmrm.Fd; 00225 ppbmrm.hist_wdist[0]=wdist; 00226 00227 /* Verbose output */ 00228 00229 if (verbose) 00230 SG_SPRINT("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n", 00231 ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, R, K); 00232 00233 /* main loop */ 00234 00235 while (ppbmrm.exitflag==0) 00236 { 00237 tstart=ttime.cur_time_diff(false); 00238 ppbmrm.nIter++; 00239 00240 /* Update H */ 00241 00242 if (ppbmrm.nCP>0) 00243 { 00244 A_2=get_cutting_plane(CPList_tail); 00245 cp_ptr=CPList_head; 00246 00247 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00248 { 00249 A_1=get_cutting_plane(cp_ptr); 00250 cp_ptr=cp_ptr->next; 00251 rsum=0.0; 00252 00253 for (uint32_t j=0; j<nDim; ++j) 00254 { 00255 rsum+=A_1[j]*A_2[j]; 00256 } 00257 00258 H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)]=rsum; 00259 } 00260 00261 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00262 { 00263 H[LIBBMRM_INDEX(ppbmrm.nCP, i, BufSize)]= 00264 H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)]; 00265 } 00266 } 00267 00268 rsum=0.0; 00269 A_2=get_cutting_plane(CPList_tail); 00270 00271 for (uint32_t i=0; i<nDim; ++i) 00272 rsum+=A_2[i]*A_2[i]; 00273 00274 H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]=rsum; 00275 00276 diag_H[ppbmrm.nCP]=H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]; 00277 I[ppbmrm.nCP]=1; 00278 00279 ppbmrm.nCP++; 00280 beta[ppbmrm.nCP]=0.0; // [beta; 0] 00281 00282 /* tune alpha cycle */ 00283 /* ---------------------------------------------------------------------- */ 00284 00285 flag=true; 00286 isThereGoodSolution=false; 00287 LIBBMRM_MEMCPY(beta_start, beta, ppbmrm.nCP*sizeof(float64_t)); 00288 LIBBMRM_MEMCPY(I_start, I, ppbmrm.nCP*sizeof(uint32_t)); 00289 qp_cnt=0; 00290 alpha_good=alpha; 00291 00292 if (tuneAlpha) 00293 { 00294 alpha_start=alpha; alpha=0.0; 00295 beta[ppbmrm.nCP]=0.0; 00296 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t)); 00297 I2[ppbmrm.nCP]=1; 00298 00299 /* add alpha-dependent terms to H, diag_h and b */ 00300 cp_ptr=CPList_head; 00301 00302 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00303 { 00304 rsum=0.0; 00305 A_1=get_cutting_plane(cp_ptr); 00306 cp_ptr=cp_ptr->next; 00307 00308 for (uint32_t j=0; j<nDim; ++j) 00309 rsum+=A_1[j]*prevW[j]; 00310 00311 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00312 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00313 00314 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00315 H2[LIBBMRM_INDEX(i, j, BufSize)]= 00316 H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00317 } 00318 00319 /* solve QP with current alpha */ 00320 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta, 00321 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00322 ppbmrm.qp_exitflag=qp_exitflag.exitflag; 00323 qp_cnt++; 00324 Fd_alpha0=-qp_exitflag.QP; 00325 00326 /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */ 00327 for (uint32_t i=0; i<nDim; ++i) 00328 { 00329 rsum=0.0; 00330 cp_ptr=CPList_head; 00331 00332 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00333 { 00334 A_1=get_cutting_plane(cp_ptr); 00335 cp_ptr=cp_ptr->next; 00336 rsum+=A_1[i]*beta[j]; 00337 } 00338 00339 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha); 00340 } 00341 00342 sq_norm_Wdiff=0.0; 00343 00344 for (uint32_t i=0; i<nDim; ++i) 00345 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]); 00346 00347 if (::sqrt(sq_norm_Wdiff) <= K) 00348 { 00349 flag=false; 00350 00351 if (alpha!=alpha_start) 00352 alphaChanged=true; 00353 } 00354 else 00355 { 00356 alpha=alpha_start; 00357 } 00358 00359 while(flag) 00360 { 00361 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t)); 00362 LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t)); 00363 I2[ppbmrm.nCP]=1; 00364 beta[ppbmrm.nCP]=0.0; 00365 00366 /* add alpha-dependent terms to H, diag_h and b */ 00367 cp_ptr=CPList_head; 00368 00369 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00370 { 00371 rsum=0.0; 00372 A_1=get_cutting_plane(cp_ptr); 00373 cp_ptr=cp_ptr->next; 00374 00375 for (uint32_t j=0; j<nDim; ++j) 00376 rsum+=A_1[j]*prevW[j]; 00377 00378 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00379 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00380 00381 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00382 H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00383 } 00384 00385 /* solve QP with current alpha */ 00386 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta, 00387 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00388 ppbmrm.qp_exitflag=qp_exitflag.exitflag; 00389 qp_cnt++; 00390 00391 /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */ 00392 for (uint32_t i=0; i<nDim; ++i) 00393 { 00394 rsum=0.0; 00395 cp_ptr=CPList_head; 00396 00397 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00398 { 00399 A_1=get_cutting_plane(cp_ptr); 00400 cp_ptr=cp_ptr->next; 00401 rsum+=A_1[i]*beta[j]; 00402 } 00403 00404 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha); 00405 } 00406 00407 sq_norm_Wdiff=0.0; 00408 for (uint32_t i=0; i<nDim; ++i) 00409 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]); 00410 00411 if (::sqrt(sq_norm_Wdiff) > K) 00412 { 00413 /* if there is a record of some good solution 00414 * (i.e. adjust alpha by division by 2) */ 00415 00416 if (isThereGoodSolution) 00417 { 00418 LIBBMRM_MEMCPY(beta, beta_good, ppbmrm.nCP*sizeof(float64_t)); 00419 LIBBMRM_MEMCPY(I2, I_good, ppbmrm.nCP*sizeof(uint32_t)); 00420 alpha=alpha_good; 00421 qp_exitflag=qp_exitflag_good; 00422 flag=false; 00423 } 00424 else 00425 { 00426 if (alpha == 0) 00427 { 00428 alpha=1.0; 00429 alphaChanged=true; 00430 } 00431 else 00432 { 00433 alpha*=2; 00434 alphaChanged=true; 00435 } 00436 } 00437 } 00438 else 00439 { 00440 if (alpha > 0) 00441 { 00442 /* keep good solution and try for alpha /= 2 if previous alpha was 1 */ 00443 LIBBMRM_MEMCPY(beta_good, beta, ppbmrm.nCP*sizeof(float64_t)); 00444 LIBBMRM_MEMCPY(I_good, I2, ppbmrm.nCP*sizeof(uint32_t)); 00445 alpha_good=alpha; 00446 qp_exitflag_good=qp_exitflag; 00447 isThereGoodSolution=true; 00448 00449 if (alpha!=1.0) 00450 { 00451 alpha/=2.0; 00452 alphaChanged=true; 00453 } 00454 else 00455 { 00456 alpha=0.0; 00457 alphaChanged=true; 00458 } 00459 } 00460 else 00461 { 00462 flag=false; 00463 } 00464 } 00465 } 00466 } 00467 else 00468 { 00469 alphaChanged=false; 00470 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t)); 00471 LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t)); 00472 00473 /* add alpha-dependent terms to H, diag_h and b */ 00474 cp_ptr=CPList_head; 00475 00476 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00477 { 00478 rsum=0.0; 00479 A_1=get_cutting_plane(cp_ptr); 00480 cp_ptr=cp_ptr->next; 00481 00482 for (uint32_t j=0; j<nDim; ++j) 00483 rsum+=A_1[j]*prevW[j]; 00484 00485 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00486 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00487 00488 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00489 H2[LIBBMRM_INDEX(i, j, BufSize)]= 00490 H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00491 } 00492 /* solve QP with current alpha */ 00493 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta, 00494 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00495 ppbmrm.qp_exitflag=qp_exitflag.exitflag; 00496 qp_cnt++; 00497 } 00498 00499 /* ----------------------------------------------------------------------------------------------- */ 00500 00501 /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */ 00502 ppbmrm.nzA=0; 00503 00504 for (uint32_t aaa=0; aaa<ppbmrm.nCP; ++aaa) 00505 { 00506 if (beta[aaa]>epsilon) 00507 { 00508 ++ppbmrm.nzA; 00509 ICPcounter[aaa]=0; 00510 } 00511 else 00512 { 00513 ICPcounter[aaa]+=1; 00514 } 00515 } 00516 00517 /* W update */ 00518 for (uint32_t i=0; i<nDim; ++i) 00519 { 00520 rsum=0.0; 00521 cp_ptr=CPList_head; 00522 00523 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00524 { 00525 A_1=get_cutting_plane(cp_ptr); 00526 cp_ptr=cp_ptr->next; 00527 rsum+=A_1[i]*beta[j]; 00528 } 00529 00530 W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha); 00531 } 00532 00533 /* risk and subgradient computation */ 00534 R = model->risk(subgrad, W); 00535 b[ppbmrm.nCP]=-R; 00536 add_cutting_plane(&CPList_tail, map, A, 00537 find_free_idx(map, BufSize), subgrad, nDim); 00538 00539 sq_norm_W=0.0; 00540 sq_norm_Wdiff=0.0; 00541 sq_norm_prevW=0.0; 00542 00543 for (uint32_t j=0; j<nDim; ++j) 00544 { 00545 b[ppbmrm.nCP]+=subgrad[j]*W[j]; 00546 sq_norm_W+=W[j]*W[j]; 00547 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]); 00548 sq_norm_prevW+=prevW[j]*prevW[j]; 00549 } 00550 00551 /* compute Fp and Fd */ 00552 ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff; 00553 ppbmrm.Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW; 00554 00555 /* gamma + tuneAlpha flag */ 00556 if (alphaChanged) 00557 { 00558 eps=1.0-(ppbmrm.Fd/ppbmrm.Fp); 00559 gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps)); 00560 } 00561 00562 if ((lastFp-ppbmrm.Fp) <= gamma) 00563 { 00564 tuneAlpha=true; 00565 } 00566 else 00567 { 00568 tuneAlpha=false; 00569 } 00570 00571 /* Stopping conditions - set only with nonzero alpha */ 00572 if (alpha==0.0) 00573 { 00574 if (ppbmrm.Fp-ppbmrm.Fd<=TolRel*LIBBMRM_ABS(ppbmrm.Fp)) 00575 ppbmrm.exitflag=1; 00576 00577 if (ppbmrm.Fp-ppbmrm.Fd<=TolAbs) 00578 ppbmrm.exitflag=2; 00579 } 00580 00581 if (ppbmrm.nCP>=BufSize) 00582 ppbmrm.exitflag=-1; 00583 00584 tstop=ttime.cur_time_diff(false); 00585 00586 /* compute wdist (= || W_{t+1} - W_{t} || ) */ 00587 sq_norm_Wdiff=0.0; 00588 00589 for (uint32_t i=0; i<nDim; ++i) 00590 { 00591 sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]); 00592 } 00593 00594 wdist=::sqrt(sq_norm_Wdiff); 00595 00596 /* Keep history of Fp, Fd, wdist */ 00597 ppbmrm.hist_Fp[ppbmrm.nIter]=ppbmrm.Fp; 00598 ppbmrm.hist_Fd[ppbmrm.nIter]=ppbmrm.Fd; 00599 ppbmrm.hist_wdist[ppbmrm.nIter]=wdist; 00600 00601 /* Verbose output */ 00602 if (verbose) 00603 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", 00604 ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, ppbmrm.Fp-ppbmrm.Fd, 00605 (ppbmrm.Fp-ppbmrm.Fd)/ppbmrm.Fp, R, ppbmrm.nCP, ppbmrm.nzA, wdist, alpha, 00606 qp_cnt, gamma, tuneAlpha); 00607 00608 /* Check size of Buffer */ 00609 if (ppbmrm.nCP>=BufSize) 00610 { 00611 ppbmrm.exitflag=-2; 00612 SG_SERROR("Buffer exceeded.\n"); 00613 } 00614 00615 /* keep w_t + Fp */ 00616 LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t)); 00617 lastFp=ppbmrm.Fp; 00618 00619 /* Inactive Cutting Planes (ICP) removal */ 00620 if (cleanICP) 00621 { 00622 /* find ICP */ 00623 cntICP=0; 00624 cntACP=0; 00625 cp_ptr=CPList_head; 00626 uint32_t tmp_idx=0; 00627 00628 while (cp_ptr != CPList_tail) 00629 { 00630 if (ICPcounter[tmp_idx++]>=cleanAfter) 00631 { 00632 ICPs[cntICP++]=cp_ptr->address; 00633 } 00634 else 00635 { 00636 ACPs[cntACP++]=tmp_idx-1; 00637 } 00638 00639 cp_ptr=cp_ptr->next; 00640 } 00641 00642 /* do ICP removal */ 00643 if (cntICP > 0) 00644 { 00645 nCP_new=ppbmrm.nCP-cntICP; 00646 00647 for (uint32_t i=0; i<cntICP; ++i) 00648 { 00649 tmp_idx=0; 00650 cp_ptr=CPList_head; 00651 00652 while(cp_ptr->address != ICPs[i]) 00653 { 00654 cp_ptr=cp_ptr->next; 00655 tmp_idx++; 00656 } 00657 00658 remove_cutting_plane(&CPList_head, &CPList_tail, map, ICPs[i]); 00659 00660 LIBBMRM_MEMMOVE(b+tmp_idx, b+tmp_idx+1, 00661 (ppbmrm.nCP-tmp_idx)*sizeof(float64_t)); 00662 LIBBMRM_MEMMOVE(beta+tmp_idx, beta+tmp_idx+1, 00663 (ppbmrm.nCP-tmp_idx)*sizeof(float64_t)); 00664 LIBBMRM_MEMMOVE(diag_H+tmp_idx, diag_H+tmp_idx+1, 00665 (ppbmrm.nCP-tmp_idx)*sizeof(float64_t)); 00666 LIBBMRM_MEMMOVE(I+tmp_idx, I+tmp_idx+1, 00667 (ppbmrm.nCP-tmp_idx)*sizeof(uint32_t)); 00668 LIBBMRM_MEMMOVE(ICPcounter+tmp_idx, ICPcounter+tmp_idx+1, 00669 (ppbmrm.nCP-tmp_idx)*sizeof(uint32_t)); 00670 } 00671 00672 /* H */ 00673 for (uint32_t i=0; i < nCP_new; ++i) 00674 { 00675 for (uint32_t j=0; j < nCP_new; ++j) 00676 { 00677 H_buff[LIBBMRM_INDEX(i, j, BufSize)]= 00678 H[LIBBMRM_INDEX(ACPs[i], ACPs[j], BufSize)]; 00679 } 00680 } 00681 00682 for (uint32_t i=0; i<nCP_new; ++i) 00683 for (uint32_t j=0; j<nCP_new; ++j) 00684 H[LIBBMRM_INDEX(i, j, BufSize)]=H_buff[LIBBMRM_INDEX(i, j, BufSize)]; 00685 00686 ppbmrm.nCP=nCP_new; 00687 } 00688 } 00689 } /* end of main loop */ 00690 00691 ppbmrm.hist_Fp.resize_vector(ppbmrm.nIter); 00692 ppbmrm.hist_Fd.resize_vector(ppbmrm.nIter); 00693 ppbmrm.hist_wdist.resize_vector(ppbmrm.nIter); 00694 00695 cp_ptr=CPList_head; 00696 00697 while(cp_ptr!=NULL) 00698 { 00699 cp_ptr2=cp_ptr; 00700 cp_ptr=cp_ptr->next; 00701 LIBBMRM_FREE(cp_ptr2); 00702 cp_ptr2=NULL; 00703 } 00704 00705 cp_list=NULL; 00706 00707 cleanup: 00708 00709 LIBBMRM_FREE(H); 00710 LIBBMRM_FREE(b); 00711 LIBBMRM_FREE(beta); 00712 LIBBMRM_FREE(A); 00713 LIBBMRM_FREE(subgrad); 00714 LIBBMRM_FREE(diag_H); 00715 LIBBMRM_FREE(I); 00716 LIBBMRM_FREE(ICPcounter); 00717 LIBBMRM_FREE(ICPs); 00718 LIBBMRM_FREE(ACPs); 00719 LIBBMRM_FREE(H_buff); 00720 LIBBMRM_FREE(map); 00721 LIBBMRM_FREE(prevW); 00722 LIBBMRM_FREE(wt); 00723 LIBBMRM_FREE(beta_start); 00724 LIBBMRM_FREE(beta_good); 00725 LIBBMRM_FREE(I_start); 00726 LIBBMRM_FREE(I_good); 00727 LIBBMRM_FREE(I2); 00728 LIBBMRM_FREE(b2); 00729 LIBBMRM_FREE(diag_H2); 00730 LIBBMRM_FREE(H2); 00731 00732 if (cp_list) 00733 LIBBMRM_FREE(cp_list); 00734 00735 return(ppbmrm); 00736 } 00737 }