/home/coin/SVN-release/OS-2.4.1/Couenne/src/cut/sdpcuts/Heuristics.cpp

Go to the documentation of this file.
00001 /* $Id: Heuristics.cpp 508 2011-02-15 21:52:44Z pbelotti $
00002  *
00003  * Name:    Heuristics.hpp
00004  * Author:  Andrea Qualizza
00005  * Purpose: 
00006  *
00007  * This file is licensed under the Eclipse Public License (EPL)
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 // constructor
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         //initialize heurLP solver interface
00076         heurLPimproveSi_.messageHandler()->setLogLevel(0);
00077         // maximization problem
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(); // we solve via primal simplex to permit a warm start later on with resolve();
00101 
00102 
00103 
00104         //initialize min norm LP Heur solver interface
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); //minimization
00113         //add variables in MNLPSi
00114         for(int i=0;i<n;i++) {
00115                 MNLPSi_.addCol (0, NULL, NULL,  xlb[i], xub[i], 0.0); //x_i
00116         }
00117         for(int i=0;i<t;i++) {
00118                 MNLPSi_.addCol (0, NULL, NULL,  ylb[i], yub[i], 0.0); //y_i
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); //r^+_f(i,0)
00124                 MNLPSi_.addCol (0, NULL, NULL,  0.0, MNLPSi_.getInfinity(), objcoeff); //r^-_f(i,0)
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); //r^+_f(i,j)
00134                         MNLPSi_.addCol (0, NULL, NULL,  0.0, MNLPSi_.getInfinity(), objcoeff); //r^-_f(i,j)
00135                 }
00136         }
00137 
00138         //add constraints in MNLPSi
00139         for (int k=0;k<cons;k++) {
00140                 // original constraints
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                 // constraints x_i +r^+_f(i,0) -r^-_f(i,0)= x_i
00149                 CoinPackedVector vector;
00150                 vector.clear();
00151                 vector.insert(i,1.0); //x_i
00152                 vector.insert(idx,1.0); //r+
00153                 idx++;
00154                 vector.insert(idx,-1.0); //r-
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                         // constraints appr.X_ij +r^+_f(i,j) -r^-_f(i,j)= X_ij
00161                         CoinPackedVector vector;
00162                         vector.clear();
00163                         vector.insert(idx,1.0); //r+
00164                         idx++;
00165                         vector.insert(idx,-1.0); //r-
00166                         idx++;
00167                         MNLPSi_.addRow(vector,'E',0.0,0);
00168                 }
00169         }
00170         MNLPSi_.initialSolve(); // we solve via primal simplex to permit a warm start later on with resolve();
00171 } // Heuristics()
00172 /***********************************************************************/
00173 // destructor
00174 Heuristics::~Heuristics() {
00175         delete [] heurLbRowAdded_;
00176         delete [] bestSol_;
00177         delete [] xxTSol_;
00178         delete [] MNSol_;
00179         delete [] temp_row_;
00180 } // ~Heuristics()
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         // xx^T heuristic solution
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         // matrix norm minimization heuristic solution
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 } // run()
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         // compute xx^T
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 } // xxTHeur()
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         //update original constraints
00287         for(int k=0;k<cons_;k++) {
00288                 for (int i=0;i<n_+t_;i++) //clear temp_row_
00289                         temp_row_[i] = 0.0;
00290                 for (int i=0;i<n_;i++)
00291                         temp_row_[i] = origMat_[k][i]; //coeff of x_i
00292                 for (int i=0;i<t_;i++)
00293                         temp_row_[N_+i] = origMat_[k][N_+i]; //coeff of y_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         //update constraints x_i +r^+_f(i,0) -r^-_f(i,0)= x_i
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 { //i != j
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         // compute xx^T
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 printf("Current Solution:\n");
00353 for (int i=0;i<n_;i++)
00354         printf("x_%d=%.2f ",i,x[i]);
00355 printf("\n");
00356 for (int i=0;i<n_;i++)
00357         for (int j=i;j<n_;j++)
00358                 printf("x_%d_%d=%.2f ",i,j,x[indexQ(i,j,n_)]);
00359 if (t_>0) {
00360 printf("\n");
00361 for (int i=0;i<t_;i++)
00362         printf("y_%d=%.2f ",i,y[i]);
00363 }
00364 printf("\n");
00365 int status = feasibility_check(n_,t_,cons_,MNSol_,origMat_,origRhs_,origSense_,xlb_,xub_,ylb_,yub_);
00366 double value = evaluateSolution(n_,t_,MNSol_,b_,c_,Q_,objConst_);
00367 printf("---->feas:%d  value=%.2f\n",status,value);
00368 
00369 printf("curr sol value=%.2f\n",evaluateSolution(n_,t_,const_cast <double *>(si_->getColSolution ()),b_,c_,Q_,objConst_));
00370 
00371 for (int i=0;i<N_+t_;i++)
00372         printf("%.2f ",const_cast <double *>(si_->getColSolution ())[i]);
00373 printf("\n");
00374 for (int i=0;i<n_;i++)
00375         printf("%.2f ",b_[i]);
00376 for (int i=0;i<n_;i++)
00377         for (int j=i;j<n_;j++)
00378                 printf("%.2f ",Q_[i][j]);
00379 for (int i=0;i<t_;i++)
00380         printf("%.2f ",c_[i]);
00381 printf("\n");
00382 MNLPSi_.writeLp("test", "lp",1e-5,10,5,0.0,true);
00383 
00384 exit(9);
00385 */
00386 /*
00387 int status = feasibility_check(n_,t_,cons_,MNSol_,origMat_,origRhs_,origSense_,xlb_,xub_,ylb_,yub_);
00388 double value = evaluateSolution(n_,t_,MNSol_,b_,c_,Q_,objConst_);
00389 printf("---->feas:%d  value=%.2f\n",status,value);
00390 */
00391         return MNSol_;
00392 } //MNHeur()
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"); //should never happen
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         // by changing the y_i variables via LP we can try to
00444         //      -improve the current solution heurSol if it is feasible
00445         //      -make the current solution heurSol feasible if it is infeasible and all the violated
00446         //       constraints contain a a variable y_i
00447         // if there are no y_i variables (t_=0), heurLP is useless
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 } // processSol()
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 } // improveSolution()
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 

Generated on Thu Nov 10 03:05:43 2011 by  doxygen 1.4.7