SHOGUN
v3.2.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 * Written (W) 2012 Viktor Gal 00008 * Copyright (C) 2008 Chun-Nam Yu 00009 */ 00010 00011 #include <shogun/structure/CCSOSVM.h> 00012 #include <shogun/mathematics/Mosek.h> 00013 00014 using namespace shogun; 00015 00016 CCCSOSVM::CCCSOSVM() 00017 : CLinearStructuredOutputMachine() 00018 { 00019 init(); 00020 } 00021 00022 CCCSOSVM::CCCSOSVM(CStructuredModel* model, SGVector<float64_t> w) 00023 : CLinearStructuredOutputMachine(model, model->get_labels()) 00024 { 00025 init(); 00026 00027 if (w.vlen) 00028 { 00029 set_w(w); 00030 } 00031 else 00032 { 00033 m_w.resize_vector(m_model->get_dim()); 00034 m_w.zero(); 00035 } 00036 } 00037 00038 CCCSOSVM::~CCCSOSVM() 00039 { 00040 #ifdef USE_MOSEK 00041 MSK_deleteenv(&m_msk_env); 00042 #endif 00043 } 00044 00045 int32_t CCCSOSVM::mosek_qp_optimize(float64_t** G, float64_t* delta, float64_t* alpha, int32_t k, float64_t* dual_obj, float64_t rho) 00046 { 00047 #ifdef USE_MOSEK 00048 int32_t t; 00049 index_t Q_size = k*(k+1)/2; 00050 SGVector<float64_t> c(k); 00051 MSKlidxt *aptrb; 00052 MSKlidxt *aptre; 00053 MSKidxt *asub; 00054 SGVector<float64_t> aval(k); 00055 MSKboundkeye bkc[1]; 00056 float64_t blc[1]; 00057 float64_t buc[1]; 00058 MSKboundkeye *bkx; 00059 SGVector<float64_t> blx(k); 00060 SGVector<float64_t> bux(k); 00061 MSKidxt *qsubi,*qsubj; 00062 SGVector<float64_t> qval(Q_size); 00063 00064 MSKtask_t task; 00065 MSKrescodee r; 00066 00067 aptrb = (MSKlidxt*) SG_MALLOC(MSKlidxt, k); 00068 aptre = (MSKlidxt*) SG_MALLOC(MSKlidxt, k); 00069 asub = (MSKidxt*) SG_MALLOC(MSKidxt, k); 00070 bkx = (MSKboundkeye*) SG_MALLOC(MSKboundkeye, k); 00071 qsubi = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size); 00072 qsubj = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size); 00073 00074 00075 /* DEBUG */ 00076 /* 00077 for (int32_t i=0;i<k;i++) 00078 printf("delta: %.4f\n", delta[i]); 00079 00080 printf("G:\n"); 00081 for (int32_t i=0;i<k;i++) 00082 { 00083 for (int32_t j=0;j<k;j++) 00084 printf("%.4f ", G[i][j]); 00085 printf("\n"); 00086 } 00087 fflush(stdout); 00088 */ 00089 /* DEBUG */ 00090 00091 for (int32_t i=0; i < k;i++) 00092 { 00093 c[i] = -delta[i]; 00094 aptrb[i] = i; 00095 aptre[i] = i+1; 00096 asub[i] = 0; 00097 aval[i] = 1.0; 00098 bkx[i] = MSK_BK_LO; 00099 blx[i] = 0.0; 00100 bux[i] = MSK_INFINITY; 00101 } 00102 bkc[0] = MSK_BK_FX; 00103 blc[0] = m_C; 00104 buc[0] = m_C; 00105 /* 00106 bkc[0] = MSK_BK_UP; 00107 blc[0] = -MSK_INFINITY; 00108 buc[0] = m_C; 00109 */ 00110 00111 /* create the optimization task */ 00112 r = MSK_maketask(m_msk_env, 1, k, &task); 00113 00114 if (r != MSK_RES_OK) 00115 SG_ERROR("Could not create MOSEK task: %d\n", r) 00116 00117 r = MSK_inputdata(task, 00118 1,k, 00119 1,k, 00120 c,0.0, 00121 aptrb,aptre, 00122 asub,aval, 00123 bkc,blc,buc, 00124 bkx,blx,bux); 00125 00126 if (r != MSK_RES_OK) 00127 SG_ERROR("Error setting input data: %d\n", r) 00128 00129 /* coefficients for the Gram matrix */ 00130 t = 0; 00131 for (int32_t i=0;i<k;i++) 00132 { 00133 for (int32_t j=0;j<=i;j++) 00134 { 00135 qsubi[t] = i; 00136 qsubj[t] = j; 00137 qval[t] = G[i][j]/(1+rho); 00138 t++; 00139 } 00140 } 00141 00142 r = MSK_putqobj(task, k*(k+1)/2, qsubi, qsubj, qval); 00143 if (r != MSK_RES_OK) 00144 SG_ERROR("Error MSK_putqobj: %d\n", r) 00145 00146 /* DEBUG */ 00147 /* 00148 printf("t: %ld\n", t); 00149 for (int32_t i=0;i<t;i++) { 00150 printf("qsubi: %d, qsubj: %d, qval: %.4f\n", qsubi[i], qsubj[i], qval[i]); 00151 } 00152 fflush(stdout); 00153 */ 00154 /* DEBUG */ 00155 00156 /* set relative tolerance gap (DEFAULT = 1E-8)*/ 00157 MSK_putdouparam(task, MSK_DPAR_INTPNT_TOL_REL_GAP, 1E-8); 00158 00159 r = MSK_optimize(task); 00160 00161 if (r != MSK_RES_OK) 00162 SG_ERROR("Error MSK_optimize: %d\n", r) 00163 00164 MSK_getsolutionslice(task, 00165 MSK_SOL_ITR, 00166 MSK_SOL_ITEM_XX, 00167 0, 00168 k, 00169 alpha); 00170 00171 /* output the objective value */ 00172 MSK_getprimalobj(task, MSK_SOL_ITR, dual_obj); 00173 00174 MSK_deletetask(&task); 00175 00176 /* free the memory */ 00177 SG_FREE(aptrb); 00178 SG_FREE(aptre); 00179 SG_FREE(asub); 00180 SG_FREE(bkx); 00181 SG_FREE(qsubi); 00182 SG_FREE(qsubj); 00183 00184 return r; 00185 #else 00186 return -1; 00187 #endif 00188 } 00189 00190 bool CCCSOSVM::train_machine(CFeatures* data) 00191 { 00192 if (data) 00193 set_features(data); 00194 00195 SGVector<float64_t> alpha; 00196 float64_t** G; /* Gram matrix */ 00197 DynArray<SGSparseVector<float64_t> > dXc; /* constraint matrix */ 00198 // DOC **dXc; /* constraint matrix */ 00199 SGVector<float64_t> delta; /* rhs of constraints */ 00200 SGSparseVector<float64_t> new_constraint; 00201 float64_t dual_obj=0, alphasum; 00202 int32_t iter, size_active; 00203 float64_t value; 00204 SGVector<int32_t> idle; /* for cleaning up */ 00205 float64_t margin; 00206 float64_t primal_obj; 00207 SGVector<float64_t> proximal_rhs; 00208 SGVector<float64_t> gammaG0; 00209 float64_t min_rho = 0.001; 00210 float64_t serious_counter=0; 00211 float64_t rho = 1.0; /* temporarily set it to 1 first */ 00212 00213 float64_t expected_descent, primal_obj_b=-1, reg_master_obj; 00214 int32_t null_step=1; 00215 float64_t kappa=0.1; 00216 float64_t temp_var; 00217 float64_t proximal_term, primal_lower_bound; 00218 00219 float64_t v_k; 00220 float64_t obj_difference; 00221 SGVector<float64_t> cut_error; // cut_error[i] = alpha_{k,i} at current center x_k 00222 float64_t sigma_k; 00223 float64_t m2 = 0.2; 00224 float64_t m3 = 0.9; 00225 float64_t gTd; 00226 float64_t last_sigma_k=0; 00227 00228 float64_t initial_primal_obj; 00229 int32_t suff_decrease_cond=0; 00230 float64_t decrease_proportion = 0.2; // start from 0.2 first 00231 00232 float64_t z_k_norm; 00233 float64_t last_z_k_norm=0; 00234 00235 /* warm start */ 00236 SGVector<float64_t> w_b = m_w.clone(); 00237 00238 iter = 0; 00239 size_active = 0; 00240 G = NULL; 00241 00242 new_constraint = find_cutting_plane(&margin); 00243 value = margin - new_constraint.dense_dot(1.0, m_w.vector, m_w.vlen, 0); 00244 00245 primal_obj_b = primal_obj = 0.5*m_w.dot(m_w.vector, m_w.vector, m_w.vlen)+m_C*value; 00246 primal_lower_bound = 0; 00247 expected_descent = -primal_obj_b; 00248 initial_primal_obj = primal_obj_b; 00249 00250 SG_INFO("Running CCCP inner loop solver: ") 00251 00252 while ((!suff_decrease_cond) && (expected_descent<-m_eps) && (iter<m_max_iter)) 00253 { 00254 ++iter; 00255 ++size_active; 00256 00257 SG_DEBUG("ITER %d\n", iter) 00258 SG_PRINT(".") 00259 00260 /* add constraint */ 00261 dXc.resize_array(size_active); 00262 dXc[size_active - 1] = new_constraint; 00263 // dXc[size_active - 1].add(new_constraint); 00264 /* 00265 dXc = (DOC**)realloc(dXc, sizeof(DOC*)*size_active); 00266 dXc[size_active-1] = (DOC*)malloc(sizeof(DOC)); 00267 dXc[size_active-1]->fvec = new_constraint; 00268 dXc[size_active-1]->slackid = 1; // only one common slackid (one-slack) 00269 dXc[size_active-1]->costfactor = 1.0; 00270 */ 00271 delta.resize_vector(size_active); 00272 delta[size_active-1] = margin; 00273 alpha.resize_vector(size_active); 00274 alpha[size_active-1] = 0.0; 00275 idle.resize_vector(size_active); 00276 idle[size_active-1] = 0; 00277 /* proximal point */ 00278 proximal_rhs.resize_vector(size_active); 00279 cut_error.resize_vector(size_active); 00280 // note g_i = - new_constraint 00281 cut_error[size_active-1] = m_C*(new_constraint.dense_dot(1.0, w_b.vector, w_b.vlen, 0) - new_constraint.dense_dot(1.0, m_w.vector, m_w.vlen, 0)); 00282 cut_error[size_active-1] += (primal_obj_b - 0.5*w_b.dot(w_b.vector, w_b.vector, w_b.vlen)); 00283 cut_error[size_active-1] -= (primal_obj - 0.5*m_w.dot(m_w.vector, m_w.vector, m_w.vlen)); 00284 00285 gammaG0.resize_vector(size_active); 00286 00287 /* update Gram matrix */ 00288 G = SG_REALLOC(float64_t*, G, size_active-1, size_active); 00289 G[size_active - 1] = NULL; 00290 for (index_t j=0; j < size_active;j++) 00291 { 00292 G[j] = SG_REALLOC(float64_t, G[j], size_active-1, size_active); 00293 } 00294 for (index_t j=0; j < size_active-1; j++) 00295 { 00296 G[size_active-1][j] = dXc[size_active-1].sparse_dot(dXc[j]); 00297 G[j][size_active-1] = G[size_active-1][j]; 00298 } 00299 G[size_active-1][size_active-1] = dXc[size_active-1].sparse_dot(dXc[size_active-1]); 00300 00301 /* update gammaG0 */ 00302 if (null_step==1) 00303 { 00304 gammaG0[size_active-1] = dXc[size_active-1].dense_dot(1.0, w_b.vector, w_b.vlen, 0); 00305 } 00306 else 00307 { 00308 for (index_t i = 0; i < size_active; i++) 00309 gammaG0[i] = dXc[i].dense_dot(1.0, w_b.vector, w_b.vlen, 0); 00310 } 00311 00312 /* update proximal_rhs */ 00313 for (index_t i = 0; i < size_active; i++) 00314 { 00315 switch(m_qp_type) 00316 { 00317 case MOSEK: 00318 proximal_rhs[i] = delta[i] - rho/(1+rho)*gammaG0[i]; 00319 break; 00320 case SVMLIGHT: 00321 proximal_rhs[i] = (1+rho)*delta[i] - rho*gammaG0[i]; 00322 break; 00323 default: 00324 SG_ERROR("Invalid QPType: %d\n", m_qp_type) 00325 } 00326 } 00327 00328 switch(m_qp_type) 00329 { 00330 case MOSEK: 00331 /* solve QP to update alpha */ 00332 dual_obj = 0; 00333 mosek_qp_optimize(G, proximal_rhs.vector, alpha.vector, size_active, &dual_obj, rho); 00334 break; 00335 case SVMLIGHT: 00336 /* TODO: port required functionality from the latest SVM^light into shogun 00337 * in order to be able to support this 00338 * 00339 if (size_active>1) 00340 { 00341 if (svmModel!=NULL) 00342 free_model(svmModel,0); 00343 svmModel = (MODEL*)my_malloc(sizeof(MODEL)); 00344 svm_learn_optimization(dXc,proximal_rhs,size_active,sm->sizePsi,&lparm,&kparm,NULL,svmModel,alpha); 00345 } 00346 else 00347 { 00348 ASSERT(size_active==1) 00349 alpha[0] = m_C; 00350 } 00351 */ 00352 break; 00353 default: 00354 SG_ERROR("Invalid QPType: %d\n", m_qp_type) 00355 } 00356 00357 /* DEBUG */ 00358 //printf("r: %d\n", r); fflush(stdout); 00359 //printf("dual: %.16lf\n", dual_obj); 00360 /* END DEBUG */ 00361 00362 m_w.zero(); 00363 for (index_t j = 0; j < size_active; j++) 00364 { 00365 if (alpha[j]>m_C*m_alpha_thrld) 00366 { 00367 // TODO: move this to SGVector 00368 // basically it's vector[i]= scale*sparsevector[i] 00369 for (index_t k = 0; k < dXc[j].num_feat_entries; k++) 00370 { 00371 index_t idx = dXc[j].features[k].feat_index; 00372 m_w[idx] += alpha[j]/(1+rho)*dXc[j].features[k].entry; 00373 } 00374 } 00375 } 00376 00377 if (m_qp_type == SVMLIGHT) 00378 { 00379 /* compute dual obj */ 00380 dual_obj = +0.5*(1+rho)*m_w.dot(m_w.vector, m_w.vector, m_w.vlen); 00381 for (int32_t j=0;j<size_active;j++) 00382 dual_obj -= proximal_rhs[j]/(1+rho)*alpha[j]; 00383 } 00384 00385 z_k_norm = CMath::sqrt(m_w.dot(m_w.vector, m_w.vector, m_w.vlen)); 00386 m_w.vec1_plus_scalar_times_vec2(m_w.vector, rho/(1+rho), w_b.vector, w_b.vlen); 00387 00388 /* detect if step size too small */ 00389 sigma_k = 0; 00390 alphasum = 0; 00391 for (index_t j = 0; j < size_active; j++) 00392 { 00393 sigma_k += alpha[j]*cut_error[j]; 00394 alphasum+=alpha[j]; 00395 } 00396 sigma_k/=m_C; 00397 gTd = -m_C*(new_constraint.dense_dot(1.0, m_w.vector, m_w.vlen, 0) 00398 - new_constraint.dense_dot(1.0, w_b.vector, w_b.vlen, 0)); 00399 00400 for (index_t j = 0; j < size_active; j++) 00401 SG_DEBUG("alpha[%d]: %.8g, cut_error[%d]: %.8g\n", j, alpha[j], j, cut_error[j]) 00402 SG_DEBUG("sigma_k: %.8g\n", sigma_k) 00403 SG_DEBUG("alphasum: %.8g\n", alphasum) 00404 SG_DEBUG("g^T d: %.8g\n", gTd) 00405 00406 /* update cleanup information */ 00407 for (index_t j = 0; j < size_active; j++) 00408 { 00409 if (alpha[j]<m_alpha_thrld*m_C) 00410 idle[j]++; 00411 else 00412 idle[j]=0; 00413 } 00414 00415 new_constraint = find_cutting_plane(&margin); 00416 value = margin - new_constraint.dense_dot(1.0, m_w.vector, m_w.vlen, 0); 00417 00418 /* print primal objective */ 00419 primal_obj = 0.5*m_w.dot(m_w.vector, m_w.vector, m_w.vlen)+m_C*value; 00420 00421 SG_DEBUG("ITER PRIMAL_OBJ %.4f\n", primal_obj) 00422 00423 temp_var = w_b.dot(w_b.vector, w_b.vector, w_b.vlen); 00424 proximal_term = 0.0; 00425 for (index_t i=0; i < m_model->get_dim(); i++) 00426 proximal_term += (m_w[i]-w_b[i])*(m_w[i]-w_b[i]); 00427 00428 reg_master_obj = -dual_obj+0.5*rho*temp_var/(1+rho); 00429 expected_descent = reg_master_obj - primal_obj_b; 00430 00431 v_k = (reg_master_obj - proximal_term*rho/2) - primal_obj_b; 00432 00433 primal_lower_bound = CMath::max(primal_lower_bound, reg_master_obj - 0.5*rho*(1+rho)*proximal_term); 00434 00435 SG_DEBUG("ITER REG_MASTER_OBJ: %.4f\n", reg_master_obj) 00436 SG_DEBUG("ITER EXPECTED_DESCENT: %.4f\n", expected_descent) 00437 SG_DEBUG("ITER PRIMLA_OBJ_B: %.4f\n", primal_obj_b) 00438 SG_DEBUG("ITER RHO: %.4f\n", rho) 00439 SG_DEBUG("ITER ||w-w_b||^2: %.4f\n", proximal_term) 00440 SG_DEBUG("ITER PRIMAL_LOWER_BOUND: %.4f\n", primal_lower_bound) 00441 SG_DEBUG("ITER V_K: %.4f\n", v_k) 00442 SG_DEBUG("ITER margin: %.4f\n", margin) 00443 SG_DEBUG("ITER psi*-psi: %.4f\n", value-margin) 00444 00445 obj_difference = primal_obj - primal_obj_b; 00446 00447 if (primal_obj<primal_obj_b+kappa*expected_descent) 00448 { 00449 /* extra condition to be met */ 00450 if ((gTd>m2*v_k)||(rho<min_rho+1E-8)) 00451 { 00452 SG_DEBUG("SERIOUS STEP\n") 00453 00454 /* update cut_error */ 00455 for (index_t i = 0; i < size_active; i++) 00456 { 00457 cut_error[i] -= (primal_obj_b - 0.5*w_b.dot(w_b.vector, w_b.vector, w_b.vlen)); 00458 cut_error[i] -= m_C*dXc[i].dense_dot(1.0, w_b.vector, w_b.vlen, 0); 00459 cut_error[i] += (primal_obj - 0.5*m_w.dot(m_w, m_w, m_w.vlen)); 00460 cut_error[i] += m_C*dXc[i].dense_dot(1.0, m_w.vector, m_w.vlen, 0); 00461 } 00462 primal_obj_b = primal_obj; 00463 /* copy w_b <- m_w */ 00464 for (index_t i=0; i < m_model->get_dim(); i++) 00465 { 00466 w_b[i] = m_w[i]; 00467 } 00468 null_step = 0; 00469 serious_counter++; 00470 } 00471 else 00472 { 00473 /* increase step size */ 00474 SG_DEBUG("NULL STEP: SS(ii) FAILS.\n") 00475 00476 serious_counter--; 00477 rho = CMath::max(rho/10,min_rho); 00478 } 00479 } 00480 else 00481 { /* no sufficient decrease */ 00482 serious_counter--; 00483 00484 if ((cut_error[size_active-1]>m3*last_sigma_k)&&(CMath::abs(obj_difference)>last_z_k_norm+last_sigma_k)) 00485 { 00486 SG_DEBUG("NULL STEP: NS(ii) FAILS.\n") 00487 rho = CMath::min(10*rho,m_max_rho); 00488 } 00489 else 00490 SG_DEBUG("NULL STEP\n") 00491 } 00492 /* update last_sigma_k */ 00493 last_sigma_k = sigma_k; 00494 last_z_k_norm = z_k_norm; 00495 00496 00497 /* break away from while loop if more than certain proportioal decrease in primal objective */ 00498 if (primal_obj_b/initial_primal_obj<1-decrease_proportion) 00499 suff_decrease_cond = 1; 00500 00501 /* clean up */ 00502 if (iter % m_cleanup_check == 0) 00503 { 00504 size_active = resize_cleanup(size_active, idle, alpha, delta, gammaG0, proximal_rhs, &G, dXc, cut_error); 00505 ASSERT(size_active == proximal_rhs.vlen) 00506 } 00507 } // end cutting plane while loop 00508 00509 SG_INFO(" Inner loop optimization finished.\n") 00510 00511 for (index_t j = 0; j < size_active; j++) 00512 SG_FREE(G[j]); 00513 SG_FREE(G); 00514 00515 /* copy */ 00516 for (index_t i=0; i < m_model->get_dim(); i++) 00517 m_w[i] = w_b[i]; 00518 00519 m_primal_obj = primal_obj_b; 00520 00521 return true; 00522 } 00523 00524 SGSparseVector<float64_t> CCCSOSVM::find_cutting_plane(float64_t* margin) 00525 { 00526 SGVector<float64_t> new_constraint(m_model->get_dim()); 00527 int32_t psi_size = m_model->get_dim(); 00528 00529 index_t num_samples = m_model->get_features()->get_num_vectors(); 00530 /* find cutting plane */ 00531 *margin = 0; 00532 new_constraint.zero(); 00533 for (index_t i = 0; i < num_samples; i++) 00534 { 00535 CResultSet* result = m_model->argmax(m_w, i); 00536 new_constraint.add(result->psi_truth); 00537 result->psi_pred.scale(-1.0); 00538 new_constraint.add(result->psi_pred); 00539 /* 00540 printf("%.16lf %.16lf\n", 00541 SGVector<float64_t>::dot(result->psi_truth.vector, result->psi_truth.vector, result->psi_truth.vlen), 00542 SGVector<float64_t>::dot(result->psi_pred.vector, result->psi_pred.vector, result->psi_pred.vlen)); 00543 */ 00544 *margin += result->delta; 00545 SG_UNREF(result); 00546 } 00547 /* scaling */ 00548 float64_t scale = 1/(float64_t)num_samples; 00549 new_constraint.scale(scale); 00550 *margin *= scale; 00551 00552 /* find the nnz elements in new_constraint */ 00553 index_t l = 0; 00554 for (index_t i=0; i < psi_size; i++) 00555 { 00556 if (CMath::abs(new_constraint[i])>1E-10) 00557 l++; // non-zero 00558 } 00559 /* TODO: does this really work good? 00560 l = CMath::get_num_nonzero(new_constraint.vector, new_constraint.vlen); 00561 */ 00562 /* create a sparse vector of the nnz of new_constraint */ 00563 SGSparseVector<float64_t> cut_plane(l); 00564 index_t k = 0; 00565 for (index_t i = 0; i < psi_size; i++) 00566 { 00567 if (CMath::abs(new_constraint[i])>1E-10) 00568 { 00569 cut_plane.features[k].feat_index = i; 00570 cut_plane.features[k].entry = new_constraint[i]; 00571 k++; 00572 } 00573 } 00574 00575 return cut_plane; 00576 } 00577 00578 int32_t CCCSOSVM::resize_cleanup(int32_t size_active, SGVector<int32_t>& idle, SGVector<float64_t>&alpha, 00579 SGVector<float64_t>& delta, SGVector<float64_t>& gammaG0, 00580 SGVector<float64_t>& proximal_rhs, float64_t ***ptr_G, 00581 DynArray<SGSparseVector<float64_t> >& dXc, SGVector<float64_t>& cut_error) 00582 { 00583 int32_t i, j, new_size_active; 00584 index_t k; 00585 00586 float64_t **G = *ptr_G; 00587 00588 i=0; 00589 while ((i<size_active)&&(idle[i]<m_idle_iter)) 00590 i++; 00591 00592 j=i; 00593 while((j<size_active)&&(idle[j]>=m_idle_iter)) 00594 j++; 00595 00596 while (j<size_active) 00597 { 00598 /* copying */ 00599 alpha[i] = alpha[j]; 00600 delta[i] = delta[j]; 00601 gammaG0[i] = gammaG0[j]; 00602 cut_error[i] = cut_error[j]; 00603 00604 SG_FREE(G[i]); 00605 G[i] = G[j]; 00606 G[j] = NULL; 00607 // free_example(dXc[i],0); 00608 dXc[i] = dXc[j]; 00609 // dXc[j] = NULL; 00610 00611 i++; 00612 j++; 00613 while((j<size_active)&&(idle[j]>=m_idle_iter)) 00614 j++; 00615 } 00616 for (k=i;k<size_active;k++) 00617 { 00618 if (G[k]!=NULL) SG_FREE(G[k]); 00619 // if (dXc[k].num_feat_entries > 0) SG_UNREF(dXc[k]); 00620 } 00621 new_size_active = i; 00622 alpha.resize_vector(new_size_active); 00623 delta.resize_vector(new_size_active); 00624 gammaG0.resize_vector(new_size_active); 00625 proximal_rhs.resize_vector(new_size_active); 00626 G = SG_REALLOC(float64_t*, G, size_active, new_size_active); 00627 dXc.resize_array(new_size_active); 00628 cut_error.resize_vector(new_size_active); 00629 00630 /* resize G and idle */ 00631 i=0; 00632 while ((i<size_active)&&(idle[i]<m_idle_iter)) 00633 i++; 00634 00635 j=i; 00636 while((j<size_active)&&(idle[j]>=m_idle_iter)) 00637 j++; 00638 00639 while (j<size_active) 00640 { 00641 idle[i] = idle[j]; 00642 for (k=0;k<new_size_active;k++) 00643 G[k][i] = G[k][j]; 00644 00645 i++; 00646 j++; 00647 while((j<size_active)&&(idle[j]>=m_idle_iter)) 00648 j++; 00649 } 00650 idle.resize_vector(new_size_active); 00651 for (k=0;k<new_size_active;k++) 00652 G[k] = SG_REALLOC(float64_t, G[k], size_active, new_size_active); 00653 00654 *ptr_G = G; 00655 00656 return new_size_active; 00657 } 00658 00659 void CCCSOSVM::init() 00660 { 00661 m_C = 1.0; 00662 m_eps = 1E-3; 00663 m_alpha_thrld = 1E-14; 00664 m_cleanup_check = 100; 00665 m_idle_iter = 20; 00666 m_max_iter = 1000; 00667 m_max_rho = m_C; 00668 m_primal_obj = CMath::INFTY; 00669 m_qp_type = MOSEK; 00670 00671 #ifdef USE_MOSEK 00672 MSKrescodee r = MSK_RES_OK; 00673 00674 /* create mosek environment */ 00675 #if (MSK_VERSION_MAJOR == 6) 00676 r = MSK_makeenv(&m_msk_env, NULL, NULL, NULL, NULL); 00677 #elif (MSK_VERSION_MAJOR == 7) 00678 r = MSK_makeenv(&m_msk_env, NULL); 00679 #else 00680 #error "Unsupported Mosek version" 00681 #endif 00682 00683 /* check return code */ 00684 if (r != MSK_RES_OK) 00685 SG_ERROR("Error while creating mosek env: %d\n", r) 00686 00687 /* initialize the environment */ 00688 r = MSK_initenv(m_msk_env); 00689 if (r != MSK_RES_OK) 00690 SG_ERROR("Error while initializing mosek env: %d\n", r) 00691 #endif 00692 00693 SG_ADD(&m_C, "m_C", "C", MS_NOT_AVAILABLE); 00694 SG_ADD(&m_eps, "m_eps", "Epsilon", MS_NOT_AVAILABLE); 00695 SG_ADD(&m_alpha_thrld, "m_alpha_thrld", "Alpha threshold", MS_NOT_AVAILABLE); 00696 SG_ADD(&m_cleanup_check, "m_cleanup_check", "Cleanup after given number of iterations", MS_NOT_AVAILABLE); 00697 SG_ADD(&m_idle_iter, "m_idle_iter", "Maximum number of idle iteration", MS_NOT_AVAILABLE); 00698 SG_ADD(&m_max_iter, "m_max_iter", "Maximum number of iterations", MS_NOT_AVAILABLE); 00699 SG_ADD(&m_max_rho, "m_max_rho", "Max rho", MS_NOT_AVAILABLE); 00700 SG_ADD(&m_primal_obj, "m_primal_obj", "Primal objective value", MS_NOT_AVAILABLE); 00701 SG_ADD((machine_int_t*) &m_qp_type, "m_qp_type", "QP Solver Type", MS_NOT_AVAILABLE); 00702 } 00703 00704 EMachineType CCCSOSVM::get_classifier_type() 00705 { 00706 return CT_CCSOSVM; 00707 }