/home/coin/SVN-release/OS-2.4.1/Bonmin/experimental/Bcp/BM_tm.cpp

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines Corporation 2006, 2007
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // Laszlo Ladanyi, International Business Machines Corporation
00007 // Pierre Bonami, Carnegie Mellon University
00008 
00009 #include "OsiClpSolverInterface.hpp"
00010 
00011 #include "BCP_problem_core.hpp"
00012 #include "BCP_tm.hpp"
00013 #include "BCP_lp.hpp"
00014 
00015 // #include "BonIpoptSolver.hpp"
00016 // #ifdef COIN_HAS_FILTERSQP
00017 // # include "BonFilterSolver.hpp"
00018 // #endif
00019 // #include "BonOACutGenerator2.hpp"
00020 // #include "BonEcpCuts.hpp"
00021 // #include "BonOaNlpOptim.hpp"
00022 #include "BonAmplSetup.hpp"
00023 
00024 #include "BM.hpp"
00025 
00026 //#############################################################################
00027 
00028 void
00029 BM_tm::readIpopt()
00030 {
00031     if (par.entry(BM_par::IpoptParamfile).length() != 0) {
00032         FILE* ipopt = fopen(par.entry(BM_par::IpoptParamfile).c_str(), "r");
00033         if (!ipopt) {
00034             throw BCP_fatal_error("Non-existent IpoptParamfile\n");
00035         }
00036         fseek(ipopt, 0, SEEK_END);
00037         int len = ftell(ipopt) + 1;
00038         fseek(ipopt, 0, SEEK_SET);
00039         char* ipopt_content = new char[len+1];
00040         len = fread(ipopt_content, 1, len, ipopt);
00041         ipopt_file_content.assign(ipopt_content, len);
00042         delete[] ipopt_content;
00043     }
00044 
00045     FILE* nl = fopen(par.entry(BM_par::NL_filename).c_str(), "r");
00046     if (!nl) {
00047         throw BCP_fatal_error("Non-existent NL_filename\n");
00048     }
00049     fseek(nl, 0, SEEK_END);
00050     int len = ftell(nl) + 1;
00051     fseek(nl, 0, SEEK_SET);
00052     char* nl_content = new char[len+1];
00053     len = fread(nl_content, 1, len, nl);
00054     nl_file_content.assign(nl_content, len);
00055     delete[] nl_content;
00056 }
00057 
00058 /****************************************************************************/
00059 
00060 void
00061 BM_tm::initialize_core(BCP_vec<BCP_var_core*>& vars,
00062                        BCP_vec<BCP_cut_core*>& cuts,
00063                        BCP_lp_relax*& matrix)
00064 {
00065     char* argv_[3];
00066     char** argv = argv_;
00067     argv[0] = NULL;
00068     argv[1] = strdup(par.entry(BM_par::NL_filename).c_str());
00069     argv[2] = NULL;
00070     
00071     /* Get the basic options. */
00072     Bonmin::BonminAmplSetup bonmin;
00073     bonmin.readOptionsFile(par.entry(BM_par::IpoptParamfile).c_str());
00074     bonmin.initialize(argv);    
00075     
00076     free(argv[1]);
00077 
00078     Bonmin::OsiTMINLPInterface& nlp = *bonmin.nonlinearSolver();
00079 #if defined(BM_DISREGARD_SOS)
00080     const Bonmin::TMINLP::SosInfo * sos = nlp.model()->sosConstraints();
00081     if (sos->num > 0) {
00082       printf("There are SOS constraints... disregarding them...\n");
00083     }
00084 #endif
00085     
00086     OsiSolverInterface& clp  = *bonmin.continuousSolver();
00087     
00088     const int numCols = clp.getNumCols();
00089     const int numRows = clp.getNumRows();
00090 
00091     const double* clb = clp.getColLower();
00092     const double* cub = clp.getColUpper();
00093 
00094     double* obj = new double[numCols];
00095     if (bonmin.getAlgorithm() == Bonmin::B_BB /* pure B&B */) {
00096       std::cout<<"Doing branch and bound"<<std::endl;
00097       CoinFillN(obj, numCols, 0.0);
00098       matrix = NULL;
00099     } else {
00100       std::cout<<"Doing hybrid"<<std::endl;
00101       CoinDisjointCopyN(clp.getObjCoefficients(), numCols, obj);
00102       cuts.reserve(numRows);
00103       const double* rlb = clp.getRowLower();
00104       const double* rub = clp.getRowUpper();
00105       for (int i = 0; i < numRows; ++i) {
00106         cuts.push_back(new BCP_cut_core(rlb[i], rub[i]));
00107       }
00108       matrix = new BCP_lp_relax(true /*column major ordered*/);
00109       matrix->copyOf(*clp.getMatrixByCol(), obj, clb, cub, rlb, rub);
00110     }
00111 
00112     vars.reserve(numCols);
00113     for (int i = 0; i < numCols; ++i)   {
00114       BCP_var_t type = BCP_ContinuousVar;
00115       if (clp.isBinary(i)) type = BCP_BinaryVar;
00116       if (clp.isInteger(i)) type = BCP_IntegerVar;
00117       vars.push_back(new BCP_var_core(type, obj[i], clb[i], cub[i]));
00118     }
00119     delete[] obj;
00120 
00121     /* Initialize pseudocosts */
00122     int nObj = 0;
00123     for (int i = 0; i < numCols; ++i) {
00124         if (nlp.isInteger(i)) {
00125             ++nObj;
00126         }
00127     }
00128 #if ! defined(BM_DISREGARD_SOS)
00129     const Bonmin::TMINLP::SosInfo * sos = nlp.model()->sosConstraints();
00130     nObj += sos->num;
00131 #endif
00132     pseudoCosts_.initialize(nObj);
00133 }
00134 
00135 /****************************************************************************/
00136 void
00137 BM_tm::create_root(BCP_vec<BCP_var*>& added_vars,
00138                    BCP_vec<BCP_cut*>& added_cuts,
00139                    BCP_user_data*& user_data)
00140 {
00141     BM_node* data = new BM_node;
00142     data->numNlpFailed_ = 0;
00143     user_data = data;
00144 }
00145 
00146 void
00147 BM_tm::write_AMPL_solution(const BCP_solution* sol,
00148                            bool write_file, bool write_screen)
00149 {
00150   const BCP_solution_generic* bg =
00151     dynamic_cast<const BCP_solution_generic*>(sol);
00152 
00153   /* Parse again the input file so that we have a nice and clean ampl
00154      setup */  
00155   char* argv_[3];
00156   char** argv = argv_;
00157   argv[0] = NULL;
00158   argv[1] = strdup(par.entry(BM_par::NL_filename).c_str());
00159   argv[2] = NULL;
00160   Bonmin::BonminAmplSetup bonmin;
00161   bonmin.initialize(argv);    
00162   
00163   Bonmin::OsiTMINLPInterface& nlpSolver = *bonmin.nonlinearSolver();
00164   free(argv[1]);
00165   OsiSolverInterface& clp = *bonmin.continuousSolver();
00166 
00167   const int numCols = clp.getNumCols();
00168 
00169   int i;
00170   
00171   /* Create a dense vector with the value of each variable. Round the
00172      integer vars (they were tested to be near enough to integrality) so
00173      in the printouts we won't have 9.99999991234, etc.) */
00174   double* dsol = new double[numCols];
00175   for (i = 0; i < numCols; ++i) {
00176     dsol[i] = 0.0;
00177   }
00178   const BCP_vec<BCP_var*>& vars = bg->_vars;
00179   const BCP_vec<double>& values = bg->_values;;
00180   const int size = vars.size();
00181   for (i = 0; i < size; ++i) {
00182     const int ind = vars[i]->bcpind();
00183     const double v = values[i];
00184     const BCP_var_t type = vars[i]->var_type();
00185     dsol[ind] = (type == BCP_ContinuousVar) ? v : floor(v+0.5);
00186   }
00187 
00188   if (write_screen) {
00189     /* Display the solution on the screen */
00190     printf("bonmin: feasible solution found.  Objective value: %f\n",
00191            bg->objective_value());
00192     for (i = 0; i < size; ++i) {
00193       printf("    index: %5i   value: %f\n",
00194              vars[i]->bcpind(), dsol[vars[i]->bcpind()]);
00195     }
00196     printf("\n");
00197   }
00198 
00199   if (write_file) {
00200     /* create the AMPL solfile */
00201     nlpSolver.model()->finalize_solution(Bonmin::TMINLP::SUCCESS, nlpSolver.getNumCols(), 
00202                                          dsol, bg->objective_value());
00203   }
00204   delete[] dsol;
00205 }
00206 
00207 //#############################################################################
00208 
00209 void
00210 BM_tm::process_message(BCP_buffer& buf)
00211 {
00212   int msgtag;
00213   buf.unpack(msgtag);
00214   if (msgtag == BM_PseudoCostUpdate) {
00215     receive_pseudo_cost_update(buf);
00216   }
00217 }
00218 
00219 /****************************************************************************/
00220 
00221 void
00222 BM_tm::display_final_information(const BCP_lp_statistics& lp_stat)
00223 {
00224   bool write_screen = false;
00225   BCP_tm_prob *p = getTmProblemPointer();
00226   if (p->param(BCP_tm_par::TmVerb_FinalStatistics)) {
00227     printf("TM: Running time: %.3f\n", CoinWallclockTime() - p->start_time);
00228     printf("TM: search tree size: %i   ( processed %i )   max depth: %i\n",
00229            int(p->search_tree.size()), int(p->search_tree.processed()),
00230            p->search_tree.maxdepth());
00231     lp_stat.display();
00232 
00233     if (! p->feas_sol) {
00234       printf("TM: No feasible solution is found\n");
00235     } else {
00236       printf("TM: The best solution found has value %f\n",
00237              p->feas_sol->objective_value());
00238       if (p->param(BCP_tm_par::TmVerb_BestFeasibleSolution)) {
00239         write_screen = true;
00240       }
00241     }
00242   }
00243   if (p->feas_sol) {
00244     write_AMPL_solution(p->feas_sol, true, write_screen);
00245   }
00246 }
00247 
00248 /****************************************************************************/
00249 
00250 void
00251 BM_tm::display_feasible_solution(const BCP_solution* sol)
00252 {
00253   // For now, we want to write also the AMPL solution file...(?)
00254   // This needs to be changed though
00255   write_AMPL_solution(sol, true, true);
00256 }
00257 
00258 struct BMSearchTreeCompareBest {
00259   static const std::string compName;
00260   static inline const char* name() { return "BMSearchTreeCompareBest"; }
00261   inline bool operator()(const CoinTreeSiblings* x,
00262                          const CoinTreeSiblings* y) const {
00263     double allowable_gap = 1e-8;
00264     const double xq = x->currentNode()->getQuality();
00265     const double yq = y->currentNode()->getQuality();
00266     const int xd = x->currentNode()->getDepth();
00267     const int yd = y->currentNode()->getDepth();
00268     return ((xq < yq-allowable_gap) ||
00269             (xq <= yq+allowable_gap && xd > yd));
00270   }
00271 };
00272 
00273 const std::string
00274 BMSearchTreeCompareBest::compName("BMSearchTreeCompareBest");
00275 
00276 void
00277 BM_tm::init_new_phase(int phase,
00278                       BCP_column_generation& colgen,
00279                       CoinSearchTreeBase*& candidates)
00280 {
00281   BCP_tm_prob* p = getTmProblemPointer();
00282   colgen = BCP_DoNotGenerateColumns_Fathom;
00283   switch (p->param(BCP_tm_par::TreeSearchStrategy)) {
00284   case BCP_BestFirstSearch:
00285     candidates = new CoinSearchTree<BMSearchTreeCompareBest>;
00286     printf("Creating candidate list with BMSearchTreeCompareBest\n");
00287     break;
00288   case BCP_BreadthFirstSearch:
00289     candidates = new CoinSearchTree<CoinSearchTreeCompareBreadth>;
00290     printf("Creating candidate list with CoinSearchTreeCompareBreadth\n");
00291     break;
00292   case BCP_DepthFirstSearch:
00293     candidates = new CoinSearchTree<CoinSearchTreeCompareDepth>;
00294     printf("Creating candidate list with CoinSearchTreeCompareDepth\n");
00295     break;
00296   case BCP_PreferredFirstSearch:
00297     candidates = new CoinSearchTree<CoinSearchTreeComparePreferred>;
00298     printf("Creating candidate list with CoinSearchTreeComparePreferred\n");
00299     break;
00300   }
00301 }
00302 
00303 //#############################################################################
00304 
00305 void
00306 BM_tm::receive_pseudo_cost_update(BCP_buffer& buf)
00307 {
00308   double* upTotalChange = pseudoCosts_.upTotalChange();
00309   double* downTotalChange = pseudoCosts_.downTotalChange();
00310   int* upNumber = pseudoCosts_.upNumber();
00311   int* downNumber = pseudoCosts_.downNumber();
00312   int objInd;
00313   int branch;
00314   double pcChange;
00315   while (true) {
00316     buf.unpack(objInd);
00317     if (objInd == -1) {
00318       break;
00319     }
00320     buf.unpack(branch);
00321     buf.unpack(pcChange);
00322     if (branch == 0) {
00323       downTotalChange[objInd] += pcChange;
00324       ++downNumber[objInd];
00325     } else {
00326       upTotalChange[objInd] += pcChange;
00327       ++upNumber[objInd];
00328     }
00329   }
00330 }
00331 
00332 //-----------------------------------------------------------------------------
00333 
00334 void
00335 BM_tm::pack_pseudo_costs(BCP_buffer& buf)
00336 {
00337   buf.pack(pseudoCosts_.upTotalChange(), pseudoCosts_.numberObjects());
00338   buf.pack(pseudoCosts_.upNumber(), pseudoCosts_.numberObjects());
00339   buf.pack(pseudoCosts_.downTotalChange(), pseudoCosts_.numberObjects());
00340   buf.pack(pseudoCosts_.downNumber(), pseudoCosts_.numberObjects());
00341 }

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