00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012
00013 #include <Heuristics.hpp>
00014 #include <CoinPackedVector.hpp>
00015 #include <OsiSolverInterface.hpp>
00016 #include <OsiXxxSolverInterface.hpp>
00017 #include <sdpcuts.hpp>
00018 #include <tracer.hpp>
00019 #include <misc_util.hpp>
00020
00021 #define indexQ(i,j,n) ((n) + (i) * (2*(n)-1-(i)) / 2 + (j))
00022
00023
00024 Heuristics::Heuristics( const int n,
00025 const int t,
00026 const int cons,
00027 const double objConst,
00028 const double *b,
00029 const double *c,
00030 const double **Q,
00031 const double **origMat,
00032 const double *origRhs,
00033 const char *origSense,
00034 const double *xlb,
00035 const double *xub,
00036 const double *ylb,
00037 const double *yub,
00038 const OsiSolverInterface *si,
00039 Tracer *tracer
00040 ):
00041
00042 n_ (n),
00043 t_ (t),
00044 cons_ (cons),
00045 objConst_ (objConst),
00046 b_ (b),
00047 c_ (c),
00048 Q_ (Q),
00049 origMat_ (origMat),
00050 origRhs_ (origRhs),
00051 origSense_ (origSense),
00052 xlb_ (xlb),
00053 xub_ (xub),
00054 ylb_ (ylb),
00055 yub_ (yub),
00056 si_ (si),
00057 currObj_ (-DBL_MAX),
00058 bestObj_ (-DBL_MAX),
00059 tracer_ (tracer)
00060 {
00061
00062 N_ = n*(n+3)/2;
00063
00064 bestSol_ = new double[N_+t_];
00065 xxTSol_ = new double[N_+t_];
00066 MNSol_ = new double[N_+t_];
00067 for (int i=0;i<N_+t_;i++) {
00068 bestSol_[i] = -DBL_MAX;
00069 xxTSol_ [i] = -DBL_MAX;
00070 MNSol_ [i] = -DBL_MAX;
00071 }
00072
00073
00074
00075
00076 heurLPimproveSi_.messageHandler()->setLogLevel(0);
00077
00078 heurLPimproveSi_.setObjSense(-1);
00079
00080 heurLbRowAdded_ = new bool[cons_];
00081
00082 for(int i=0;i<cons_;i++)
00083 heurLbRowAdded_[i] = false;
00084 for (int i=0;i<t;i++)
00085 heurLPimproveSi_.addCol (0, NULL, NULL, ylb[i], yub[i], c[i]);
00086 for (int k=0;k<cons;k++) {
00087 double rhs = origRhs[k];
00088 CoinPackedVector vector;
00089 vector.clear();
00090 for (int j=0;j<t;j++) {
00091 double val = origMat[k][N_+j];
00092 if (val != 0.0)
00093 vector.insert(j,val);
00094 }
00095 if (vector.getNumElements() > 0) {
00096 heurLPimproveSi_.addRow(vector,origSense[k],rhs,0);
00097 heurLbRowAdded_[k] = true;
00098 }
00099 }
00100 heurLPimproveSi_.initialSolve();
00101
00102
00103
00104
00105
00106 double *coeff = const_cast <double *> ( si_->getObjCoefficients() );
00107
00108 temp_row_ = new double[N_+t];
00109 for(int i=0;i<N_+t;i++)
00110 temp_row_[i] = 0.0;
00111 MNLPSi_.messageHandler()->setLogLevel(0);
00112 MNLPSi_.setObjSense(1);
00113
00114 for(int i=0;i<n;i++) {
00115 MNLPSi_.addCol (0, NULL, NULL, xlb[i], xub[i], 0.0);
00116 }
00117 for(int i=0;i<t;i++) {
00118 MNLPSi_.addCol (0, NULL, NULL, ylb[i], yub[i], 0.0);
00119 }
00120 for(int i=0;i<n;i++) {
00121 double objcoeff;
00122 objcoeff = fabs(coeff[i]);
00123 MNLPSi_.addCol (0, NULL, NULL, 0.0, MNLPSi_.getInfinity(), objcoeff);
00124 MNLPSi_.addCol (0, NULL, NULL, 0.0, MNLPSi_.getInfinity(), objcoeff);
00125 }
00126 for (int i=0;i<n;i++) {
00127 for (int j=i;j<n;j++) {
00128 double objcoeff;
00129 if (i==j)
00130 objcoeff = fabs(coeff[indexQ(i,i,n_)]);
00131 else
00132 objcoeff = fabs(coeff[indexQ(i,j,n_)]);
00133 MNLPSi_.addCol (0, NULL, NULL, 0.0, MNLPSi_.getInfinity(), objcoeff);
00134 MNLPSi_.addCol (0, NULL, NULL, 0.0, MNLPSi_.getInfinity(), objcoeff);
00135 }
00136 }
00137
00138
00139 for (int k=0;k<cons;k++) {
00140
00141 CoinPackedVector vector;
00142 vector.clear();
00143 MNLPSi_.addRow(vector,origSense[k],origRhs[k],0);
00144 }
00145 int idx;
00146 idx = n_+t_;
00147 for (int i=0;i<n;i++) {
00148
00149 CoinPackedVector vector;
00150 vector.clear();
00151 vector.insert(i,1.0);
00152 vector.insert(idx,1.0);
00153 idx++;
00154 vector.insert(idx,-1.0);
00155 idx++;
00156 MNLPSi_.addRow(vector,'E',0.0,0);
00157 }
00158 for (int i=0;i<n;i++) {
00159 for (int j=i;j<n;j++) {
00160
00161 CoinPackedVector vector;
00162 vector.clear();
00163 vector.insert(idx,1.0);
00164 idx++;
00165 vector.insert(idx,-1.0);
00166 idx++;
00167 MNLPSi_.addRow(vector,'E',0.0,0);
00168 }
00169 }
00170 MNLPSi_.initialSolve();
00171 }
00172
00173
00174 Heuristics::~Heuristics() {
00175 delete [] heurLbRowAdded_;
00176 delete [] bestSol_;
00177 delete [] xxTSol_;
00178 delete [] MNSol_;
00179 delete [] temp_row_;
00180 }
00181
00182 int Heuristics::run() {
00183 Timer heur_timer;
00184 heur_timer.start();
00185
00186 int status;
00187 bool feasible_solution = false;
00188 double *tmpSol;
00189
00190 currObj_ = -DBL_MAX;
00191
00192
00193
00194
00195 #ifdef HEUR_XXT
00196 Timer xxtheur_timer;
00197 xxtheur_timer.start();
00198 double xxtheurvalue = -DBL_MAX;
00199 double xxt_heurlp_value = -DBL_MAX;
00200
00201 tmpSol = xxTHeur();
00202 #ifdef HEUR_IMPROVE_SOLUTION
00203 status = processSol(tmpSol,true,&xxtheurvalue,&xxt_heurlp_value);
00204 #else
00205 status = processSol(tmpSol,false,&xxtheurvalue,&xxt_heurlp_value);
00206 #endif
00207 if (status == 0)
00208 feasible_solution = true;
00209
00210 tracer_->setHeuristicsxxTTime(xxtheur_timer.time());
00211 if (xxtheurvalue > -DBL_MAX) {
00212 tracer_->setHeuristicsxxTSolution(xxtheurvalue);
00213 }
00214 if (xxt_heurlp_value > -DBL_MAX) {
00215 tracer_->setHeuristicsxxTSolutionLPHeuristicImprovement
00216 (fabs(xxt_heurlp_value-xxtheurvalue));
00217 }
00218 #endif //HEUR_XXT
00219
00220
00221
00222
00223 #ifdef HEUR_MIN_MATRIX_NORM
00224 Timer MNheur_timer;
00225 MNheur_timer.start();
00226 double MNheurvalue = -DBL_MAX;
00227 double MN_heurlp_value = -DBL_MAX;
00228
00229 tmpSol = MNHeur();
00230 #ifdef HEUR_IMPROVE_SOLUTION
00231 status = processSol(tmpSol,true,&MNheurvalue,&MN_heurlp_value);
00232 #else
00233 status = processSol(tmpSol,false,&MNheurvalue,&MN_heurlp_value);
00234 #endif
00235 if (status == 0)
00236 feasible_solution = true;
00237
00238 tracer_->setHeuristicsMNLPTime(MNheur_timer.time());
00239 if (MNheurvalue > -DBL_MAX)
00240 tracer_->setHeuristicsMNLPSolution(MNheurvalue);
00241 if (MN_heurlp_value > -DBL_MAX)
00242 tracer_->setHeuristicsMNLPSolutionLPHeuristicImprovement
00243 (fabs(MN_heurlp_value-MNheurvalue));
00244 #endif //HEUR_MIN_MATRIX_NORM
00245
00246 tracer_->setHeuristicsTime(heur_timer.time());
00247 if (bestObj_ > -DBL_MAX) {
00248 tracer_->setHeuristicsBestSolution(bestObj_);
00249 }
00250 if (feasible_solution) {
00251 tracer_->setHeuristicsCurrentSolution(currObj_);
00252 return 0;
00253 }
00254 else {
00255 return 1;
00256 }
00257 }
00258
00259 double* Heuristics::xxTHeur() {
00260
00261 double *sol = const_cast <double *> ( si_->getColSolution () );
00262 double *x,*y;
00263 x = sol;
00264 if (t_ > 0)
00265 y = &sol[N_];
00266
00267
00268 for(int i=0;i<n_;i++)
00269 xxTSol_[i] = x[i];
00270 for(int i=0;i<t_;i++)
00271 xxTSol_[i+N_] = y[i];
00272 for(int i=0;i<n_;i++)
00273 for(int j=i;j<n_;j++)
00274 xxTSol_[indexQ(i,j,n_)] = x[i] * x[j];
00275 return xxTSol_;
00276 }
00277
00278 double* Heuristics::MNHeur() {
00279 double *sol = const_cast <double *> ( si_->getColSolution () );
00280
00281 double *x,*y;
00282 x = sol;
00283 if (t_ > 0)
00284 y = &sol[N_];
00285
00286
00287 for(int k=0;k<cons_;k++) {
00288 for (int i=0;i<n_+t_;i++)
00289 temp_row_[i] = 0.0;
00290 for (int i=0;i<n_;i++)
00291 temp_row_[i] = origMat_[k][i];
00292 for (int i=0;i<t_;i++)
00293 temp_row_[N_+i] = origMat_[k][N_+i];
00294 for (int i=0;i<n_;i++)
00295 for (int j=i;j<n_;j++) {
00296 double value = sqrt(fabs(x[indexQ(i,j,n_)]))/2.0;
00297 temp_row_[i] += value;
00298 temp_row_[j] += value;
00299 }
00300
00301 for(int i=0;i<N_+t_;i++)
00302 MNLPSi_.XxxModifyCoefficient(k,i,temp_row_[i]);
00303 }
00304
00305
00306 int rowidx;
00307 rowidx = cons_;
00308 for (int i=0;i<n_;i++) {
00309 MNLPSi_.setRowBounds(rowidx+i,x[i],x[i]);
00310 }
00311 rowidx = cons_ + n_;
00312 for (int i=0;i<n_;i++)
00313 for (int j=i;j<n_;j++) {
00314 if (i==j) {
00315 double value = sqrt(fabs(x[indexQ(i,j,n_)]));
00316 MNLPSi_.XxxModifyCoefficient(rowidx,i,value);
00317 } else {
00318 double value = sqrt(fabs(x[indexQ(i,j,n_)]))/2.0;
00319 MNLPSi_.XxxModifyCoefficient(rowidx,i,value);
00320 MNLPSi_.XxxModifyCoefficient(rowidx,j,value);
00321 }
00322 MNLPSi_.setRowBounds(rowidx,x[indexQ(i,j,n_)],x[indexQ(i,j,n_)]);
00323 rowidx++;
00324 }
00325
00326 MNLPSi_.resolve();
00327 if (MNLPSi_.isAbandoned()) {
00328 printf("Heur:: HeurLP: Numerical difficulties in Solver\n");
00329 return NULL;
00330 }
00331 if (MNLPSi_.isProvenPrimalInfeasible()) {evaluateSolution(n_,t_,MNSol_,b_,c_,Q_,objConst_);
00332 #if (defined TRACE_ALL)||(defined TRACE_HEUR)
00333 printf("Heur:: MNLPHeur Problem is infeasible\n");
00334 #endif
00335 return NULL;
00336 }
00337
00338 const double *lpheurSol = MNLPSi_.getColSolution();
00339
00340
00341 for(int i=0;i<n_;i++)
00342 MNSol_[i] = lpheurSol[i];
00343 for(int i=0;i<t_;i++)
00344 MNSol_[N_+i] = lpheurSol[n_+i];
00345 for(int i=0;i<n_;i++)
00346 for(int j=i;j<n_;j++)
00347 MNSol_[indexQ(i,j,n_)] = MNSol_[i] * MNSol_[j];
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 return MNSol_;
00392 }
00393
00394 int Heuristics::processSol(double* sol, bool improveHeurLP, double *origvalue, double *lpheurvalue) {
00395 int status;
00396 double value = -DBL_MAX;
00397
00398 if (sol == NULL)
00399 return -1;
00400 status = feasibility_check(n_,t_,cons_,sol,origMat_,origRhs_,origSense_,xlb_,xub_,ylb_,yub_);
00401 if (status == FEAS_CHECK_NO_VIOLATION)
00402 value = evaluateSolution(n_,t_,sol,b_,c_,Q_,objConst_);
00403 (*origvalue) = value;
00404 #if (defined TRACE_ALL)||(defined TRACE_HEUR)
00405 if (value > -DBL_MAX)
00406 printf("Heur:: Heuristic solution value = %.3f",value);
00407 else
00408 printf("Heur:: Heuristic solution infeasible. ");
00409 switch(status) {
00410 case FEAS_CHECK_NO_VIOLATION:
00411 printf("(feasible)\n");
00412 break;
00413 case FEAS_CHECK_CONSTRAINT_VIOLATION:
00414 printf("Enforcing feasibility...\n");
00415 break;
00416 case FEAS_CHECK_CONSTRAINT_VIOLATION_NO_RECOVER:
00417 printf("Cannot enforce feasibility!\n");
00418 break;
00419 case FEAS_CHECK_BOUNDS_VIOLATION:
00420 printf("\nERROR: variable bounds are violated!!!\n");
00421 break;
00422 default:
00423 printf("\nERROR: feasibility_check() return code unknown!!\n");
00424 }
00425 #endif
00426 #ifdef TRACE_ALL
00427 printf("Heur:: Heuristic Solution:\nx = ");
00428 for(int i=0;i<n_;i++)
00429 printf("%.2f ",sol[i]);
00430 if (t_>0) {
00431 printf("\ny = ");
00432 for(int i=0;i<t_;i++)
00433 printf("%.2f ",sol[i+N_]);
00434 printf("\n");
00435 }
00436 #endif
00437 if ( (status != FEAS_CHECK_NO_VIOLATION) && (status != FEAS_CHECK_CONSTRAINT_VIOLATION) )
00438 return status;
00439
00440 if (status == FEAS_CHECK_NO_VIOLATION)
00441 update(sol,value);
00442
00443
00444
00445
00446
00447
00448 if ((improveHeurLP) && (t_>0) && (heurLP_improveSolution(sol) == 0)) {
00449 status = feasibility_check(n_,t_,cons_,sol,origMat_,origRhs_,origSense_,xlb_,xub_,ylb_,yub_);
00450 double newValue = evaluateSolution(n_,t_,sol,b_,c_,Q_,objConst_);
00451 if (status == FEAS_CHECK_NO_VIOLATION) {
00452 (*lpheurvalue) = newValue;
00453 update(sol,newValue);
00454 #if (defined TRACE_ALL)||(defined TRACE_HEUR)
00455 if (newValue > value)
00456 printf("Heur:: HeurLP improved feasible solution found! value=%.2f\n",newValue);
00457 else
00458 printf("Heur:: HeurLP could not improve the current heuristic solution\n");
00459 #endif
00460 }
00461 }
00462 return status;
00463
00464 }
00465
00466 int Heuristics::heurLP_improveSolution(double* sol) {
00467 int rowcnt = 0;
00468 for (int k=0;k<cons_;k++) {
00469 if (heurLbRowAdded_[k]) {
00470 double rhs = origRhs_[k];
00471 for (int i=0;i<N_;i++) {
00472 rhs -= origMat_[k][i] * sol[i];
00473 }
00474 if (origSense_[k] == 'G')
00475 heurLPimproveSi_.setRowLower(rowcnt,rhs);
00476 if (origSense_[k] == 'L')
00477 heurLPimproveSi_.setRowUpper(rowcnt,rhs);
00478 if (origSense_[k] == 'E') {
00479 heurLPimproveSi_.setRowLower(rowcnt,rhs);
00480 heurLPimproveSi_.setRowUpper(rowcnt,rhs);
00481 }
00482 rowcnt++;
00483 }
00484 }
00485
00486 heurLPimproveSi_.resolve();
00487
00488 if (heurLPimproveSi_.isAbandoned()) {
00489 printf("Heur:: HeurLP: Numerical difficulties in Solver\n");
00490 return -1;
00491 }
00492 if (heurLPimproveSi_.isProvenPrimalInfeasible()) {
00493 #if (defined TRACE_ALL)||(defined TRACE_HEUR)
00494 printf("Heur:: HeurLP Problem is infeasible\n");
00495 #endif
00496 return 1;
00497 }
00498
00499 const double *lpheurSol = heurLPimproveSi_.getColSolution();
00500 for(int i=0;i<t_;i++)
00501 sol[i+N_] = lpheurSol[i];
00502
00503 return 0;
00504 }
00505
00506 int Heuristics::update(double* sol, double value) {
00507 if (currObj_ < value)
00508 currObj_ = value;
00509 if (bestObj_ < value) {
00510 bestObj_ = value;
00511 for(int i=0;i<N_+t_;i++)
00512 bestSol_[i] = sol[i];
00513 return 0;
00514 } else
00515 return 1;
00516 }
00517
00518
00519
00520