/home/coin/SVN-release/OS-2.4.2/Couenne/src/interfaces/CouenneMINLPInterface.cpp

Go to the documentation of this file.
00001 /* $Id: CouenneMINLPInterface.cpp 530 2011-03-12 16:52:06Z pbelotti $
00002  *
00003  * Name:    CouenneMINLPInterface
00004  * Author:  Pietro Belotti
00005  * Purpose: MINLP interface with computation of gradient and Jacobian
00006  *          through CouenneExpressions
00007  *
00008  * (C) Pietro Belotti 2010
00009  * This file is licensed under the Eclipse Public License (EPL)
00010  */
00011 
00012 #include "CouenneMINLPInterface.hpp"
00013 
00014 using namespace Couenne;
00015 
00017 void CouenneMINLPInterface::setInitSol (const CouNumber *sol) {
00018   
00019 }
00020 
00023 CouNumber CouenneMINLPInterface::solve (CouNumber *solution) {
00024 
00025   return 0.;
00026 }
00027 
00028 
00029 #if 0
00030 //#include "BonminConfig.h"
00031 
00032 //#include "BonOsiTMINLPInterface.hpp"
00033 #include "CoinTime.hpp"
00034 #include <climits>
00035 #include <string>
00036 #include <sstream>
00037 #include "BonAuxInfos.hpp"
00038 
00039 #include "Ipopt/BonIpoptSolver.hpp"
00040 #include "Ipopt/BonIpoptWarmStart.hpp"
00041 
00042 #ifdef COIN_HAS_FILTERSQP
00043 #include "Filter/BonFilterSolver.hpp"
00044 #include "Filter/BonFilterWarmStart.hpp"
00045 //#include "Filter/BonBqpdWarmStart.hpp"
00046 #endif
00047 
00048 #include "OsiBranchingObject.hpp"
00049 #include "BonStrongBranchingSolver.hpp"
00050 
00051 //Macros to try and make messages definition less heavy
00052 #include "BonMsgUtils.hpp"
00053 
00054 using namespace Ipopt;
00055 using namespace Couenne;
00056 
00058 static void
00059 register_general_options
00060 (SmartPtr<RegisteredOptions> roptions)
00061 {
00062   roptions->SetRegisteringCategory("nlp interface option", RegisteredOptions::BonminCategory);
00063   roptions->AddStringOption3("nlp_solver",
00064                              "Choice of the solver for local optima of continuous nlp's",
00065                              "Ipopt",
00066                              "Ipopt", "Interior Point OPTimizer (https://projects.coin-or.org/Ipopt)",
00067                              "filterSQP", "Sequential quadratic programming trust region "
00068                              "algorithm (http://www-unix.mcs.anl.gov/~leyffer/solvers.html)",
00069                              "all", "run all available solvers at each node",
00070                              "Note that option will work only if the specified solver has been installed. Ipopt will usualy be installed with Bonmin by default. For FilterSQP please see http://www-unix.mcs.anl.gov/~leyffer/solvers.html on how to obtain it and https://projects.coin-or.org/Bonmin/wiki/HintTricks on how to configure Bonmin to use it.");
00071   roptions->setOptionExtraInfo("nlp_solver",15);
00072   roptions->AddBoundedIntegerOption("nlp_log_level",
00073                                     "specify NLP solver interface log level (independent from ipopt print_level).",
00074                                     0,2,1,
00075                                     "Set the level of output of the CouenneMINLPInterface : "
00076                                     "0 - none, 1 - normal, 2 - verbose"
00077                                     );
00078   roptions->setOptionExtraInfo("nlp_log_level",15);
00079 
00080   roptions->AddStringOption2("file_solution",
00081                              "Write a file bonmin.sol with the solution",
00082                              "no",
00083                              "yes","",
00084                              "no","","");
00085 
00086   roptions->AddStringOption3("warm_start",
00087                              "Select the warm start method",
00088                              "none",
00089                              "none","No warm start",
00090                              "optimum","Warm start with direct parent optimum",
00091                              "interior_point","Warm start with an interior point of direct parent",
00092                              "This will affect the function getWarmStart(), and as a consequence the warm starting in the various algorithms.");
00093   roptions->setOptionExtraInfo("warm_start",8);
00094 
00095   roptions->SetRegisteringCategory("Nlp solution robustness", RegisteredOptions::BonminCategory);
00096 
00097   roptions->AddLowerBoundedNumberOption("max_random_point_radius",
00098                                         "Set max value r for coordinate of a random point.",
00099                                         0.,1,1e5,
00100                                         "When picking a random point coordinate i will be in the interval [min(max(l,-r),u-r), max(min(u,r),l+r)] "
00101                                         "(where l is the lower bound for the variable and u is its upper bound)");
00102   roptions->setOptionExtraInfo("max_random_point_radius",8);
00103 
00104   roptions->AddStringOption3("random_point_type","method to choose a random starting point",
00105                              "Jon",
00106                              "Jon", "Choose random point uniformly between the bounds",
00107                              "Andreas", "perturb the starting point of the problem within a prescribed interval",
00108                              "Claudia", "perturb the starting point using the perturbation radius suffix information",
00109                              "");
00110   roptions->setOptionExtraInfo("random_point_type",8);
00111 
00112   roptions->AddLowerBoundedNumberOption("random_point_perturbation_interval",
00113                                         "Amount by which starting point is perturbed when choosing to pick random point by perturbating starting point",
00114                                         0.,true, 1.,
00115                                         "");
00116   roptions->setOptionExtraInfo("random_point_perturbation_interval",8);
00117                                            
00118 
00119   roptions->AddLowerBoundedIntegerOption
00120     ("num_iterations_suspect",
00121      "Number of iterations over which a node is considered \"suspect\" (for debugging purposes only, see detailed documentation).",
00122      -1,-1,
00123      "When the number of iterations to solve a node is above this number, the subproblem at this"
00124      " node is considered to be suspect and it will be outputed in a file (set to -1 to deactivate this).");
00125   roptions->setOptionExtraInfo("num_iterations_suspect",15);
00126 
00127   
00128 
00129   roptions->AddLowerBoundedIntegerOption("num_retry_unsolved_random_point",
00130                                          "Number $k$ of times that the algorithm will try to resolve an unsolved NLP with a random starting point "
00131                                          "(we call unsolved an NLP for which Ipopt is not "
00132                                          "able to guarantee optimality within the specified tolerances).",
00133                                          0,0,
00134                                          "When Ipopt fails to solve a continuous NLP sub-problem, if $k > 0$, the algorithm will "
00135                                          "try again to solve the failed NLP with $k$ new randomly chosen starting points "
00136                                          " or until the problem is solved with success.");
00137   roptions->setOptionExtraInfo("num_retry_unsolved_random_point",15);
00138 
00139 
00140   roptions->SetRegisteringCategory("Options for non-convex problems", RegisteredOptions::BonminCategory);
00141 
00142 
00143   roptions->AddLowerBoundedIntegerOption("num_resolve_at_root",
00144                                          "Number $k$ of tries to resolve the root node with different starting points.",
00145                                          0,0,
00146                                          "The algorithm will solve the root node with $k$ random starting points"
00147                                          " and will keep the best local optimum found.");
00148   roptions->setOptionExtraInfo("num_resolve_at_root",8);
00149 
00150   roptions->AddLowerBoundedIntegerOption("num_resolve_at_node",
00151                                          "Number $k$ of tries to resolve a node (other than the root) of the tree with different starting point.",
00152                                          0,0,
00153                                          "The algorithm will solve all the nodes with $k$ different random starting points "
00154                                          "and will keep the best local optimum found.");
00155   roptions->setOptionExtraInfo("num_resolve_at_node",8);
00156 
00157   roptions->AddLowerBoundedIntegerOption("num_resolve_at_infeasibles",
00158                                          "Number $k$ of tries to resolve an infeasible node (other than the root) of the tree with different starting point.",
00159                                          0,0,
00160                                          "The algorithm will solve all the infeasible nodes with $k$ different random starting points "
00161                                          "and will keep the best local optimum found.");
00162   roptions->setOptionExtraInfo("num_resolve_at_infeasibles",8);
00163 
00164 
00165 
00166 }
00167 
00168 static void register_OA_options
00169 (SmartPtr<RegisteredOptions> roptions)
00170 {
00171   roptions->SetRegisteringCategory("Outer Approximation cuts generation", RegisteredOptions::BonminCategory);
00172   
00173   roptions->AddStringOption2("disjunctive_cut_type",
00174                              "Determine if and what kind of disjunctive cuts should be computed.",
00175                              "none",
00176                              "none", "No disjunctive cuts.",
00177                              "most-fractional", "If discrete variables present, compute disjunction for most-fractional variable");
00178   roptions->setOptionExtraInfo("disjunctive_cut_type",7);
00179 
00180   roptions->AddStringOption2("oa_cuts_scope","Specify if OA cuts added are to be set globally or locally valid",
00181                              "global",
00182                              "local","Cuts are treated as globally valid",
00183                              "global", "Cuts are treated as locally valid",
00184                              "");
00185   roptions->setOptionExtraInfo("oa_cuts_scope",7);
00186 
00187   roptions->AddStringOption2("add_only_violated_oa","Do we add all OA cuts or only the ones violated by current point?",
00188                              "no",
00189                              "no","Add all cuts",
00190                              "yes","Add only violated Cuts","");
00191   roptions->setOptionExtraInfo("add_only_violated_oa",7);
00192 
00193   roptions->AddStringOption4("cut_strengthening_type",
00194                              "Determines if and what kind of cut strengthening should be performed.",
00195                              "none",
00196                              "none", "No strengthening of cuts.",
00197                              "sglobal", "Strengthen global cuts.",
00198                              "uglobal-slocal", "Unstrengthened global and strengthened local cuts",
00199                              "sglobal-slocal", "Strengthened global and strengthened local cuts",
00200                              "");
00201   roptions->setOptionExtraInfo("cut_strengthening_type",7);
00202   
00203   roptions->AddLowerBoundedNumberOption("tiny_element","Value for tiny element in OA cut",
00204                                         -0.,0,1e-08,
00205                                         "We will remove \"cleanly\" (by relaxing cut) an element lower"
00206                                         " than this.");
00207   roptions->setOptionExtraInfo("tiny_element",7);
00208 
00209   roptions->AddLowerBoundedNumberOption("very_tiny_element","Value for very tiny element in OA cut",
00210                                         -0.,0,1e-17,
00211                                         "Algorithm will take the risk of neglecting an element lower"
00212                                         " than this.");
00213   roptions->setOptionExtraInfo("very_tiny_element",7);
00214 
00215   roptions->AddLowerBoundedIntegerOption("oa_cuts_log_level",
00216                                          "level of log when generating OA cuts.",
00217                                          0, 0,
00218                                          "0: outputs nothing,\n"
00219                                          "1: when a cut is generated, its violation and index of row from which it originates,\n"
00220                                          "2: always output violation of the cut.\n"
00221                                          "3: output generated cuts incidence vectors.");
00222   roptions->setOptionExtraInfo("oa_cuts_log_level",7);
00223 
00224 }
00225 
00226 
00228 void
00229 CouenneMINLPInterface::registerOptions
00230 (SmartPtr<RegisteredOptions> roptions)
00231 {
00232   // We try to register the options - if those have been registered
00233   // already, we catch the exception and don't need to do it again
00234   try {
00235     register_general_options(roptions);
00236     register_OA_options(roptions);
00237 #ifdef COIN_HAS_FILTERSQP
00238     FilterSolver::RegisterOptions(roptions);
00239 #endif
00240 #ifdef COIN_HAS_IPOPT
00241     IpoptSolver::RegisterOptions(roptions);
00242 #endif
00243   }   
00244   catch(RegisteredOptions::OPTION_ALREADY_REGISTERED) {
00245     // skipping
00246   }
00247 
00248 
00249 }
00250 
00252 void 
00253 CouenneMINLPInterface::setAppDefaultOptions(SmartPtr<OptionsList> Options)
00254 {}
00255 
00256 CouenneMINLPInterface::Messages::Messages
00257 ():CoinMessages((int)OSITMINLPINTERFACE_DUMMY_END)
00258 {
00259   strcpy(source_ ,"NLP");
00260   ADD_MSG(SOLUTION_FOUND, std_m, 2,
00261           "After %d tries found a solution of %g (previous best %g).");
00262 
00263   ADD_MSG(INFEASIBLE_SOLUTION_FOUND, std_m, 2,
00264           "After %d tries found an solution of %g infeasible problem.");
00265 
00266   ADD_MSG(UNSOLVED_PROBLEM_FOUND, std_m, 2,
00267           "After %d tries found an solution of %g unsolved problem.");
00268   ADD_MSG(WARN_SUCCESS_WS, warn_m, 2,
00269           "Problem not solved with warm start but solved without");
00270 
00271   ADD_MSG(WARNING_RESOLVING, warn_m,2,
00272           "Trying to resolve NLP with different starting point (%d attempts).");
00273   ADD_MSG(WARN_SUCCESS_RANDOM, warn_m, 1,
00274           "Problem initially not solved but solved with a random starting point (success on %d attempt)");
00275   ADD_MSG(WARN_CONTINUING_ON_FAILURE, warn_m, 1,
00276           "Warning : continuing branching, while there are unrecovered failures at nodes");
00277 
00278   ADD_MSG(SUSPECT_PROBLEM, std_m, 2, "NLP number %d is suspect (see bounds and start file)");
00279   ADD_MSG(IPOPT_SUMMARY, std_m, 2, "Ipopt return (for %s): status %2d, iter count %4d, time %g");
00280   ADD_MSG(BETTER_SOL, std_m, 2, "Solution of value %g found on %d'th attempt");
00281 
00282   ADD_MSG(LOG_HEAD, std_m, 1,
00283           "\n          "
00284           "    Num      Status      Obj             It       time");
00285   ADD_MSG(LOG_FIRST_LINE, std_m, 1, 
00286           "    %-8d %-11s %-23.16g %-8d %-3g");
00287   ADD_MSG(LOG_LINE, std_m, 1, " %c  r%-7d %-11s %-23.16g %-8d %-3g");
00288 
00289   ADD_MSG(ALTERNATE_OBJECTIVE, std_m, 1, "Objective value recomputed with alternate objective: %g.");
00290   
00291   ADD_MSG(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, warn_m, 1, 
00292           "resolve called before any call to initialSol  can not use warm starts.");
00293   ADD_MSG(ERROR_NO_TNLPSOLVER, warn_m, 1,"Can not parse options when no IpApplication has been created");
00294   ADD_MSG(WARNING_NON_CONVEX_OA, warn_m, 1,
00295           "OA on non-convex constraint is very experimental.");                          
00296   ADD_MSG(SOLVER_DISAGREE_STATUS, warn_m, 1, "%s says problem %s, %s says %s.");
00297   ADD_MSG(SOLVER_DISAGREE_VALUE, warn_m, 1, "%s gives objective %.16g, %s gives %.16g.");
00298 
00299 }
00300 
00301 
00302 void  
00303 CouenneMINLPInterface::OaMessageHandler::print(OsiRowCut &row){
00304   FILE * fp = filePointer();
00305   const int &n = row.row().getNumElements();
00306   fprintf(fp,"Row cut has %d elements. Lower bound: %g, upper bound %g.\n", n, row.lb(), row.ub());
00307   const int * idx = row.row().getIndices();
00308   const double * val = row.row().getElements();
00309   for(int i = 0 ; i < n ; i++){
00310     fprintf(fp,"%g, x%d",val[i], idx[i]);
00311     if(i && i % 7 == 0)
00312       fprintf(fp,"\n");
00313   } 
00314 }
00315 
00316 CouenneMINLPInterface::OaMessages::OaMessages(): CoinMessages((int) OA_MESSAGES_DUMMY_END){
00317   strcpy(source_,"OaCg");
00318   ADD_MSG(VIOLATED_OA_CUT_GENERATED, std_m, 1,"Row %d, cut violation is %g: Outer approximation cut generated.");
00319 
00320   ADD_MSG(CUT_NOT_VIOLATED_ENOUGH, std_m, 2,"Row %d, cut violation is %g: Outer approximation cut not generated.");
00321 
00322   ADD_MSG(OA_CUT_GENERATED, std_m, 1,"Row %d: Outer approximation cut not generated.");
00323 }
00324 bool CouenneMINLPInterface::hasPrintedOptions=0;
00325 
00327 // Constructors and desctructors                                  //
00330 CouenneMINLPInterface::CouenneMINLPInterface():
00331   OsiSolverInterface(),
00332   tminlp_(NULL),
00333   problem_(NULL),
00334   problem_to_optimize_(NULL),
00335   feasibility_mode_(false),
00336   app_(NULL),
00337   debug_apps_(),
00338   testOthers_(false),
00339   warmstart_(NULL),
00340   rowsense_(NULL),
00341   rhs_(NULL),
00342   rowrange_(NULL),
00343   reducedCosts_(NULL),
00344   OsiDualObjectiveLimit_(1e200),
00345   hasVarNamesFile_(true),
00346   nCallOptimizeTNLP_(0),
00347   totalNlpSolveTime_(0),
00348   totalIterations_(0),
00349   maxRandomRadius_(1e08),
00350   randomGenerationType_(0),
00351   max_perturbation_(COIN_DBL_MAX),
00352   pushValue_(1e-02),
00353   numRetryInitial_(-1),
00354   numRetryResolve_(-1),
00355   numRetryInfeasibles_(-1),
00356   numRetryUnsolved_(1),
00357   messages_(),
00358   pretendFailIsInfeasible_(0),
00359   hasContinuedAfterNlpFailure_(false),
00360   numIterationSuspect_(-1),
00361   hasBeenOptimized_(false),
00362   obj_(NULL),
00363   feasibilityProblem_(NULL),
00364   jRow_(NULL),
00365   jCol_(NULL),
00366   jValues_(NULL),
00367   nnz_jac(0),
00368   constTypes_(NULL),
00369   nNonLinear_(0),
00370   tiny_(1e-8),
00371   veryTiny_(1e-20),
00372   infty_(1e100),
00373   exposeWarmStart_(false),
00374   firstSolve_(true),
00375   cutStrengthener_(NULL),
00376   oaMessages_(),
00377   oaHandler_(NULL)
00378 {
00379   oaHandler_ = new OaMessageHandler;
00380   oaHandler_->setLogLevel(0);
00381 }
00382 
00383 void 
00384 CouenneMINLPInterface::initialize(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
00385                                   Ipopt::SmartPtr<Ipopt::OptionsList> options,
00386                                   Ipopt::SmartPtr<Ipopt::Journalist> journalist,
00387                                   Ipopt::SmartPtr<TMINLP> tminlp){
00388   if(!IsValid(app_))
00389     createApplication(roptions, options, journalist);
00390   setModel(tminlp); 
00391 }
00392 
00393 void CouenneMINLPInterface::setSolver(Ipopt::SmartPtr<TNLPSolver> app){
00394   app_ = app;
00395 }
00396 
00397 
00398 void
00399 CouenneMINLPInterface::createApplication(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
00400                                          Ipopt::SmartPtr<Ipopt::OptionsList> options,
00401                                          Ipopt::SmartPtr<Ipopt::Journalist> journalist
00402                                          )
00403 {
00404   assert(!IsValid(app_));
00405   int ival;
00406   options->GetEnumValue("nlp_solver", ival, "bonmin.");
00407   Solver s = (Solver) ival;
00408   if(s == EFilterSQP){
00409     testOthers_ = false;;
00410 #ifdef COIN_HAS_FILTERSQP
00411     app_ = new Bonmin::FilterSolver(roptions, options, journalist);
00412 #else
00413     throw SimpleError("createApplication",
00414                       "Bonmin not configured to run with FilterSQP.");
00415 #endif    
00416   }
00417   else if(s == EIpopt){
00418     testOthers_ = false;
00419 #ifdef COIN_HAS_IPOPT
00420     app_ = new IpoptSolver(roptions, options, journalist);
00421 #else
00422     throw SimpleError("createApplication",
00423                       "Bonmin not configured to run with Ipopt.");
00424 #endif
00425   }
00426   else if(s == EAll){
00427 #ifdef COIN_HAS_FILTERSQP
00428     app_ = new Bonmin::FilterSolver(roptions, options, journalist);
00429 #else
00430     throw SimpleError("createApplication",
00431                       "Bonmin not configured to run with Ipopt.");
00432 #endif
00433 #ifdef COIN_HAS_IPOPT
00434     debug_apps_.push_back(new IpoptSolver(roptions, options, journalist)); 
00435 #endif
00436     testOthers_ = true;
00437   }
00438   if (!app_->Initialize("")) {
00439     throw CoinError("Error during initialization of app_","createApplication", "CouenneMINLPInterface");
00440   }
00441   for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin() ; 
00442       i != debug_apps_.end() ; i++){
00443     (*i)->Initialize("");
00444   }
00445   extractInterfaceParams();
00446   
00447 }
00448 
00450 void
00451 CouenneMINLPInterface::setModel(SmartPtr<TMINLP> tminlp)
00452 {
00453   assert(IsValid(tminlp));
00454   tminlp_ = tminlp;
00455   problem_ = new TMINLP2TNLP(tminlp_);
00456   feasibilityProblem_ = new TNLP2FPNLP
00457     (SmartPtr<TNLP>(GetRawPtr(problem_)));
00458   if(feasibility_mode_){
00459     problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
00460   }
00461   else {
00462     problem_to_optimize_ = GetRawPtr(problem_);
00463   }
00464 }
00465 
00466 
00467 
00468 void
00469 CouenneMINLPInterface::readOptionFile(const std::string & fileName)
00470 {
00471   if(IsValid(app_)){
00472     std::ifstream is;
00473     if (fileName != "") {
00474       try {
00475         is.open(fileName.c_str());
00476       }
00477       catch(std::bad_alloc) {
00478         throw CoinError("Not enough memory to open option file.\n", "readOptionFile", "CouenneMINLPInterface");
00479       }
00480     }
00481     options()->ReadFromStream(*app_->journalist(), is);
00482     extractInterfaceParams();
00483   }
00484 }
00485 
00487 CouenneMINLPInterface::CouenneMINLPInterface (const CouenneMINLPInterface &source):
00488   OsiSolverInterface(source),
00489   tminlp_(source.tminlp_),
00490   problem_(NULL),
00491   problem_to_optimize_(NULL),
00492   feasibility_mode_(source.feasibility_mode_),
00493   rowsense_(NULL),
00494   rhs_(NULL),
00495   rowrange_(NULL),
00496   reducedCosts_(NULL),
00497   OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
00498   hasVarNamesFile_(source.hasVarNamesFile_),
00499   nCallOptimizeTNLP_(0),
00500   totalNlpSolveTime_(0),
00501   totalIterations_(0),
00502   maxRandomRadius_(source.maxRandomRadius_),
00503   randomGenerationType_(source.randomGenerationType_),
00504   max_perturbation_(source.max_perturbation_),
00505   pushValue_(source.pushValue_),
00506   numRetryInitial_(source.numRetryInitial_),
00507   numRetryResolve_(source.numRetryResolve_),
00508   numRetryInfeasibles_(source.numRetryInfeasibles_),
00509   numRetryUnsolved_(source.numRetryUnsolved_),
00510   messages_(),
00511   pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
00512   hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
00513   numIterationSuspect_(source.numIterationSuspect_),
00514   hasBeenOptimized_(source.hasBeenOptimized_),
00515   obj_(NULL),
00516   feasibilityProblem_(NULL),
00517   jRow_(NULL),
00518   jCol_(NULL),
00519   jValues_(NULL),
00520   nnz_jac(source.nnz_jac),
00521   constTypes_(NULL),
00522   //    constTypesNum_(NULL),
00523   nNonLinear_(0),
00524   tiny_(source.tiny_),
00525   veryTiny_(source.veryTiny_),
00526   infty_(source.infty_),
00527   exposeWarmStart_(source.exposeWarmStart_),
00528   firstSolve_(true),
00529   cutStrengthener_(source.cutStrengthener_),
00530   oaMessages_(),
00531   oaHandler_(NULL),
00532   strong_branching_solver_(source.strong_branching_solver_)
00533 {
00534   if(defaultHandler()){
00535     messageHandler()->setLogLevel(source.messageHandler()->logLevel());
00536   }
00537   //Pass in message handler
00538   if(source.messageHandler())
00539     passInMessageHandler(source.messageHandler());
00540   // Copy options from old application
00541   if(IsValid(source.tminlp_)) {
00542     problem_ = source.problem_->clone();
00543     feasibilityProblem_ = new TNLP2FPNLP
00544       (SmartPtr<TNLP>(GetRawPtr(problem_)), source.feasibilityProblem_);
00545     if(feasibility_mode_)
00546       problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
00547     else
00548       problem_to_optimize_ = GetRawPtr(problem_);
00549     pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
00550 
00551     setAuxiliaryInfo(source.getAuxiliaryInfo());
00552     // Copy options from old application
00553     app_ = source.app_->clone();
00554     for(std::list<Ipopt::SmartPtr<TNLPSolver> >::const_iterator i = source.debug_apps_.begin();
00555         i != source.debug_apps_.end() ; i++){
00556       debug_apps_.push_back((*i)->clone());
00557     }
00558     testOthers_ = source.testOthers_;
00559   }
00560   else {
00561     throw SimpleError("Don't know how to copy an empty IpoptInterface.",
00562                       "copy constructor");
00563   }
00564 
00565   warmstart_ = source.warmstart_ ? source.warmstart_->clone() : NULL;
00566 
00567   if(source.obj_) {
00568     obj_ = new double[source.getNumCols()];
00569     CoinCopyN(source.obj_, source.getNumCols(), obj_);
00570   }
00571 
00572 
00573   oaHandler_ = new OaMessageHandler(*source.oaHandler_);;
00574 
00575 }
00576 
00577 OsiSolverInterface * 
00578 CouenneMINLPInterface::clone(bool copyData ) const
00579 {
00580   if(copyData)
00581     return new CouenneMINLPInterface(*this);
00582   else return new CouenneMINLPInterface;
00583 }
00584 
00586 CouenneMINLPInterface & CouenneMINLPInterface::operator=(const CouenneMINLPInterface& rhs)
00587 {
00588   if(this!= &rhs) {
00589     OsiSolverInterface::operator=(rhs);
00590     OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
00591     nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
00592     totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
00593     totalIterations_ = rhs.totalIterations_;
00594     maxRandomRadius_ = rhs.maxRandomRadius_;
00595     hasVarNamesFile_ = rhs.hasVarNamesFile_;
00596     pushValue_ = rhs.pushValue_;
00597 
00598     delete warmstart_;
00599     warmstart_ = NULL;
00600 
00601     if(IsValid(rhs.tminlp_)) {
00602 
00603       tminlp_ = rhs.tminlp_;
00604       problem_ = new TMINLP2TNLP(tminlp_);
00605       problem_to_optimize_ = GetRawPtr(problem_);
00606       app_ = rhs.app_->clone();
00607 
00608       warmstart_ = rhs.warmstart_ ? rhs.warmstart_->clone() : NULL;
00609 
00610       feasibilityProblem_ = new TNLP2FPNLP
00611         (SmartPtr<TNLP>(GetRawPtr(problem_)));
00612       nnz_jac = rhs.nnz_jac;
00613 
00614       if(constTypes_ != NULL) {
00615         delete [] constTypes_;
00616         constTypes_ = NULL;
00617       }
00618       if(rhs.constTypes_ != NULL) {
00619         constTypes_ = new TNLP::LinearityType[getNumRows()];
00620         CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
00621       }
00622       /*
00623         if(constTypesNum_ != NULL) {
00624         delete [] constTypesNum_;
00625         constTypesNum_ = NULL;
00626         }
00627         if(rhs.constTypesNum_ != NULL) {
00628         constTypesNum_ = new int[getNumRows()];
00629         CoinCopyN(rhs.constTypesNum_, getNumRows(), constTypesNum_);
00630         }
00631       */
00632       if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
00633         jValues_ = new double [nnz_jac];
00634         jCol_    = new Index [nnz_jac];
00635         jRow_    = new Index [nnz_jac];
00636         CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
00637         CoinCopyN(rhs.jCol_    , nnz_jac,jCol_    );
00638         CoinCopyN(rhs.jRow_    , nnz_jac,jRow_    );
00639       }
00640       else if(nnz_jac > 0) {
00641         throw CoinError("Arrays for storing jacobian are inconsistant.",
00642                         "copy constructor",
00643                         "IpoptOAInterface");
00644       }
00645       tiny_ = rhs.tiny_;
00646       veryTiny_ = rhs.veryTiny_;
00647       infty_ = rhs.infty_;
00648       exposeWarmStart_ = rhs.exposeWarmStart_;
00649 
00650     }
00651     else {
00652       tminlp_ =NULL;
00653       problem_ = NULL;
00654       app_ = NULL;
00655       feasibilityProblem_ = NULL;
00656     }
00657 
00658 
00659     if(obj_) {
00660       delete [] obj_;
00661       obj_ = NULL;
00662     }
00663     if(rhs.obj_) {
00664       obj_ = new double[rhs.getNumCols()];
00665       CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
00666     }
00667 
00668     hasVarNamesFile_ = rhs.hasVarNamesFile_;
00669 
00670     nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
00671     totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
00672     totalIterations_ = rhs.totalIterations_;
00673     maxRandomRadius_ = rhs.maxRandomRadius_;
00674     pushValue_ = rhs.pushValue_;
00675     numRetryInitial_ = rhs.numRetryInitial_;
00676     numRetryResolve_ = rhs.numRetryResolve_;
00677     numRetryInfeasibles_ = rhs.numRetryInfeasibles_;
00678     numRetryUnsolved_ = rhs.numRetryUnsolved_;
00679     pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
00680     numIterationSuspect_ = rhs.numIterationSuspect_;
00681 
00682     hasBeenOptimized_ = rhs.hasBeenOptimized_;
00683     cutStrengthener_ = rhs.cutStrengthener_;
00684 
00685     delete oaHandler_;
00686     oaHandler_ = new OaMessageHandler(*rhs.oaHandler_);
00687     strong_branching_solver_ = rhs.strong_branching_solver_;
00688 
00689     freeCachedData();
00690   }
00691   return *this;
00692 }
00693 
00694 const SmartPtr<OptionsList> CouenneMINLPInterface::options() const
00695 {
00696   if(!IsValid(app_)) {
00697     messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
00698     return NULL;
00699   }
00700   else
00701     return app_->options();
00702 }
00703 
00704 SmartPtr<OptionsList> CouenneMINLPInterface::options()
00705 {
00706   if(!IsValid(app_)) {
00707     messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
00708     return NULL;
00709   }
00710   else
00711     return app_->options();
00712 }
00713 
00715 CouenneMINLPInterface::~CouenneMINLPInterface ()
00716 {
00717   freeCachedData();
00718   delete [] jRow_;
00719   delete [] jCol_;
00720   delete [] jValues_;
00721   delete [] constTypes_;
00722   delete [] obj_;
00723   delete oaHandler_;
00724   delete warmstart_;
00725 }
00726 
00727 void
00728 CouenneMINLPInterface::freeCachedColRim()
00729 {
00730   if(reducedCosts_!=NULL) {
00731     delete []  reducedCosts_;
00732     reducedCosts_ = NULL;
00733   }
00734 
00735 }
00736 
00737 void
00738 CouenneMINLPInterface::freeCachedRowRim()
00739 {
00740   if(rowsense_!=NULL) {
00741     delete []  rowsense_;
00742     rowsense_ = NULL;
00743   }
00744   if(rhs_!=NULL) {
00745     delete []  rhs_;
00746     rhs_ = NULL;
00747   }
00748   if(rowrange_!=NULL) {
00749     delete []  rowrange_;
00750     rowrange_ = NULL;
00751   }
00752   //   if(dualsol_!=NULL)
00753   //     {
00754   //       delete []  dualsol_; dualsol_ = NULL;
00755   //     }
00756 }
00757 
00758 void
00759 CouenneMINLPInterface::freeCachedData()
00760 {
00761   freeCachedColRim();
00762   freeCachedRowRim();
00763 }
00764 
00765 const char * CouenneMINLPInterface::OPT_SYMB="OPT";
00766 const char * CouenneMINLPInterface::FAILED_SYMB="FAILED";
00767 const char * CouenneMINLPInterface::UNBOUND_SYMB="UNBOUNDED";
00768 const char * CouenneMINLPInterface::INFEAS_SYMB="INFEAS";
00770 // WarmStart Information                                                                           //
00772 
00773 
00774 void
00775 CouenneMINLPInterface::resolveForCost(int numsolve, bool keepWarmStart)
00776 {
00777   // This method assumes that a problem has just been solved and we try for a
00778   // different solution. So disregard (in fact, clear out) any warmstart info
00779   // we might have, and acquire a new one before returning.
00780   delete warmstart_;
00781   warmstart_ = NULL;
00782  
00783   Coin::SmartPtr<SimpleReferencedPtr<CoinWarmStart> > ws_backup = NULL;
00784   if(!exposeWarmStart_ && keepWarmStart){
00785     //if warm start is not exposed, we need to store the current starting point to
00786     //restore it at the end of the method
00787     ws_backup = make_referenced(app_->getUsedWarmStart(problem_));
00788   }
00790   vector<double> point(getNumCols()*3+ getNumRows());
00791   double bestBound = getObjValue();
00792   CoinCopyN(getColSolution(),
00793             getNumCols(), point());
00794   CoinCopyN(getRowPrice(),
00795             2*getNumCols()+ getNumRows(),
00796             point() + getNumCols());
00797 
00798   if(isProvenOptimal())
00799     messageHandler()->message(SOLUTION_FOUND,
00800                               messages_)
00801       <<1<<getObjValue()<<bestBound
00802       <<CoinMessageEol;
00803   else
00804     messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
00805                               messages_)
00806       <<1
00807       <<CoinMessageEol;
00808   for(int f = 0; f < numsolve ; f++) {
00809     messageHandler()->message(WARNING_RESOLVING,
00810                               messages_)
00811       <<f+1<< CoinMessageEol ;
00812     randomStartingPoint();
00813     solveAndCheckErrors(0,0,"resolveForCost");
00814 
00815 
00816     char c=' ';
00817     //Is solution better than previous
00818     if(isProvenOptimal() &&
00819        getObjValue()<bestBound) {
00820       c='*';
00821       messageHandler()->message(BETTER_SOL, messages_)<<getObjValue()<<f+1<< CoinMessageEol;
00822       CoinCopyN(getColSolution(),
00823                 getNumCols(), point());
00824       CoinCopyN(getRowPrice(),
00825                 2*getNumCols()+ getNumRows(),
00826                 point() + getNumCols());
00827       bestBound = getObjValue();
00828     }
00829 
00830     messageHandler()->message(LOG_LINE, messages_)
00831       <<c<<f+1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
00832 
00833 
00834     if(isProvenOptimal())
00835       messageHandler()->message(SOLUTION_FOUND,
00836                                 messages_)
00837         <<f+2<<getObjValue()<<bestBound
00838         <<CoinMessageEol;
00839     else if(!isAbandoned())
00840       messageHandler()->message(UNSOLVED_PROBLEM_FOUND,
00841                                 messages_)
00842         <<f+2
00843         <<CoinMessageEol;
00844     else
00845       messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
00846                                 messages_)
00847         <<f+2
00848         <<CoinMessageEol;
00849   }
00850   setColSolution(point());
00851   setRowPrice(point() + getNumCols());
00852   app_->enableWarmStart();
00853 
00854   optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
00855   hasBeenOptimized_ = true;
00856 
00857   if(!exposeWarmStart_ && keepWarmStart) {
00858     app_->setWarmStart(ws_backup->ptr(), problem_);
00859   }
00860 }
00861 
00862 void
00863 CouenneMINLPInterface::resolveForRobustness(int numsolve)
00864 {
00865   // This method assumes that a problem has just been solved and we try for a
00866   // different solution. So disregard (in fact, clear out) any warmstart info
00867   // we might have, and acquire a new one before returning.
00868   delete warmstart_;
00869   warmstart_ = NULL;
00870   
00871 
00872   CoinWarmStart * ws_backup = NULL;
00873   if(!exposeWarmStart_){
00874     //if warm start is not exposed, we need to store the current starting point to
00875     //restore it at the end of the method
00876     ws_backup = app_->getUsedWarmStart(problem_);
00877   }
00878   //std::cerr<<"Resolving the problem for robustness"<<std::endl;
00879   //First remove warm start point and resolve
00880   app_->disableWarmStart();
00881   problem()->resetStartingPoint();
00882   messageHandler()->message(WARNING_RESOLVING,
00883                             messages_)
00884     <<1<< CoinMessageEol ;
00885   solveAndCheckErrors(0,0,"resolveForRobustness");
00886 
00887 
00888   char c='*';
00889   if(isAbandoned()) {
00890     c=' ';
00891   }
00892   messageHandler()->message(LOG_LINE, messages_)
00893     <<c<<1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<
00894     app_->CPUTime()<<CoinMessageEol;
00895 
00896 
00897   if(!isAbandoned()) {
00898     messageHandler()->message(WARN_SUCCESS_WS, messages_) << CoinMessageEol ;
00899     // re-enable warmstart and get it
00900     app_->enableWarmStart();
00901     if (!exposeWarmStart_) {
00902       app_->setWarmStart(ws_backup, problem_);
00903       delete ws_backup;
00904     }
00905     return; //we won go on
00906   }
00907 
00908   //still unsolved try again with different random starting points
00909   for(int f = 0; f < numsolve ; f++) {
00910     messageHandler()->message(WARNING_RESOLVING,
00911                               messages_)
00912       <<f+2<< CoinMessageEol ;
00913 
00914     randomStartingPoint();
00915     solveAndCheckErrors(0,0,"resolveForRobustness");
00916 
00917 
00918     messageHandler()->message(IPOPT_SUMMARY, messages_)
00919       <<"resolveForRobustness"<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
00920 
00921 
00922     char c='*';
00923     if(isAbandoned()) {
00924       c=' ';
00925     }
00926     messageHandler()->message(LOG_LINE, messages_)
00927       <<c<<f+2<<statusAsString()<<getObjValue()
00928       <<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
00929 
00930 
00931     if(!isAbandoned()) {
00932       messageHandler()->message(WARN_SUCCESS_RANDOM, messages_)
00933         << f+2 << CoinMessageEol ;
00934       // re-enable warmstart and get it
00935       app_->enableWarmStart();
00936       if (!exposeWarmStart_) {
00937         app_->setWarmStart(ws_backup, problem_);
00938         delete ws_backup;
00939       }
00940         
00941       return; //we have found a solution and continue
00942     }
00943   }
00944 
00945 
00946   if(!exposeWarmStart_){
00947     app_->setWarmStart(ws_backup, problem_);
00948     delete ws_backup;
00949   }
00950   if(pretendFailIsInfeasible_) {
00951     if(pretendFailIsInfeasible_ == 1) {
00952       messageHandler()->message(WARN_CONTINUING_ON_FAILURE,
00953                                 messages_)
00954         <<CoinMessageEol;
00955       hasContinuedAfterNlpFailure_ = 1;
00956     }
00957     return;
00958   }
00959   else {
00960     std::string probName;
00961     getStrParam(OsiProbName,probName);
00962     throw newUnsolvedError(app_->errorCode(), problem_,
00963                            probName);
00964   }
00965   // Do NOT get warmstart in other cases
00966 }
00967 
00969 // Problem information methods                                    //
00972 int CouenneMINLPInterface::getNumCols() const
00973 {
00974 
00975   return problem_->num_variables();
00976 }
00977 
00978 
00980 int
00981 CouenneMINLPInterface::getNumRows() const
00982 {
00983   return problem_->num_constraints();
00984 }
00985 
00986 const double *
00987 CouenneMINLPInterface::getColLower() const
00988 {
00989   return problem_->x_l();
00990 }
00991 
00992 const double *
00993 CouenneMINLPInterface::getColUpper() const
00994 {
00995   return problem_->x_u();
00996 }
00997 
00998 #if 1
00999 
01000 
01002 const OsiSolverInterface::OsiNameVec& 
01003 CouenneMINLPInterface::getVarNames() {
01004   return getColNames();
01005 }
01006 #endif
01007 
01008 
01009 void CouenneMINLPInterface::extractSenseRhsAndRange() const
01010 {
01011   assert(rowsense_==NULL&&rhs_==NULL&&rowrange_==NULL);
01012   int numrows = problem_->num_constraints();
01013   if(numrows == 0) return;
01014   const double * rowLower = getRowLower();
01015   const double * rowUpper = getRowUpper();
01016   rowsense_ = new char [numrows];
01017   rhs_ = new double [numrows];
01018   rowrange_ = new double [numrows];
01019   for(int i = 0 ; i < numrows ; i++) {
01020     rowrange_[i]=0.;
01021     convertBoundToSense(rowLower[i], rowUpper[i], rowsense_[i], rhs_[i], rowrange_[i]);
01022   }
01023 }
01024 
01034 const char *
01035 CouenneMINLPInterface::getRowSense() const
01036 {
01037   if(rowsense_==NULL) {
01038     extractSenseRhsAndRange();
01039   }
01040   return rowsense_;
01041 }
01042 
01051 const double *
01052 CouenneMINLPInterface::getRightHandSide() const
01053 {
01054   if(rhs_==NULL) {
01055     extractSenseRhsAndRange();
01056   }
01057   return rhs_;
01058 }
01059 
01068 const double *
01069 CouenneMINLPInterface::getRowRange() const
01070 {
01071   if(rowrange_==NULL) {
01072     extractSenseRhsAndRange();
01073   }
01074   return rowrange_;
01075 }
01076 
01077 const double *
01078 CouenneMINLPInterface::getRowLower() const
01079 {
01080   return problem_->g_l();
01081 }
01082 
01083 const double *
01084 CouenneMINLPInterface::getRowUpper() const
01085 {
01086   return problem_->g_u();
01087 }
01088 
01090 bool
01091 CouenneMINLPInterface::isContinuous(int colNumber) const
01092 {
01093   return (problem_->var_types()[colNumber]==TMINLP::CONTINUOUS);
01094 }
01095 
01097 bool
01098 CouenneMINLPInterface::isBinary(int colNumber) const
01099 {
01100   return (problem_->var_types()[colNumber]==TMINLP::BINARY);
01101 }
01102 
01107 bool
01108 CouenneMINLPInterface::isInteger(int colNumber) const
01109 {
01110   return ((problem_->var_types()[colNumber]==TMINLP::BINARY)||
01111           (problem_->var_types()[colNumber]==TMINLP::INTEGER));
01112 }
01113 
01115 bool
01116 CouenneMINLPInterface::isIntegerNonBinary(int colNumber) const
01117 {
01118   return (problem_->var_types()[colNumber]==TMINLP::INTEGER);
01119 }
01121 bool
01122 CouenneMINLPInterface::isFreeBinary(int colNumber) const
01123 {
01124   return ((problem_->var_types()[colNumber]==TMINLP::BINARY)
01125           &&((getColUpper()[colNumber]-getColLower()[colNumber]) > 1 - 1e-09));
01126 }
01127 
01129 double
01130 CouenneMINLPInterface::getInfinity() const
01131 {
01132   return COIN_DBL_MAX;
01133 }
01134 
01136 const double *
01137 CouenneMINLPInterface::getColSolution() const
01138 {
01139   if(hasBeenOptimized_)
01140     return problem_->x_sol();
01141   else
01142     return problem_->x_init();
01143 }
01144 
01146 const double *
01147 CouenneMINLPInterface::getRowPrice() const
01148 {
01149   if(hasBeenOptimized_)
01150     return problem_->duals_sol();
01151   else
01152     return problem_->duals_init();
01153 }
01154 
01156 const double *
01157 CouenneMINLPInterface::getReducedCost() const
01158 {
01159   (*handler_)<<"WARNING : trying to access reduced cost in Ipopt always retrun 0"<<CoinMessageEol;
01160   if(reducedCosts_==NULL) {
01161     reducedCosts_ = new double [getNumCols()];
01162     CoinFillN(reducedCosts_,getNumCols(),0.);
01163   }
01164   return reducedCosts_;
01165 }
01166 
01169 const double *
01170 CouenneMINLPInterface::getRowActivity() const
01171 {
01172   return problem_->g_sol();
01173 }
01174 
01178 int
01179 CouenneMINLPInterface::getIterationCount() const
01180 {
01181   return app_->IterationCount();
01182 }
01183 
01184 
01187 void
01188 CouenneMINLPInterface::setColLower( int elementIndex, double elementValue )
01189 {
01190   //  if(fabs(problem_->x_l()[elementIndex]-elementValue)>1e-06)
01191   problem_->SetVariableLowerBound(elementIndex,elementValue);
01192   hasBeenOptimized_ = false;
01193 }
01194 
01197 void
01198 CouenneMINLPInterface::setColUpper( int elementIndex, double elementValue )
01199 {
01200   //  if(fabs(problem_->x_u()[elementIndex]-elementValue)>1e-06)
01201   problem_->SetVariableUpperBound(elementIndex,elementValue);
01202   hasBeenOptimized_ = false;
01203 }
01204 
01207 void
01208 CouenneMINLPInterface::setColLower( const double* array )
01209 {
01210   problem_->SetVariablesLowerBounds(problem_->num_variables(),
01211                                     array);
01212   hasBeenOptimized_ = false;
01213 }
01214 
01217 void
01218 CouenneMINLPInterface::setColUpper( const double* array )
01219 {
01220   problem_->SetVariablesUpperBounds(problem_->num_variables(), 
01221                                     array);
01222   hasBeenOptimized_ = false;
01223 }
01224 
01227 void
01228 CouenneMINLPInterface::setRowLower( int elementIndex, double elementValue )
01229 {
01230   throw SimpleError("Not implemented yet but should be if necessary.",
01231                     "setRowLower");
01232   hasBeenOptimized_ = false;
01233 }
01234 
01237 void
01238 CouenneMINLPInterface::setRowUpper( int elementIndex, double elementValue )
01239 {
01240   throw SimpleError("Not implemented yet but should be if necessary.",
01241                     "setRowUpper");
01242   hasBeenOptimized_ = false;
01243 }
01244 
01246 void
01247 CouenneMINLPInterface::setRowType(int index, char sense, double rightHandSide,
01248                                   double range)
01249 {
01250   throw SimpleError("Not implemented yet but should be if necessary.",
01251                     "setRowType");
01252   hasBeenOptimized_ = false;
01253 }
01254 
01255 
01258 void
01259 CouenneMINLPInterface::setObjSense(double s)
01260 {
01261   throw SimpleError("Can not change objective sense of an Ipopt problem.",
01262                     "setObjSense");
01263   hasBeenOptimized_ = false;
01264 }
01265 
01275 void
01276 CouenneMINLPInterface::setColSolution(const double *colsol)
01277 {
01278   problem_->setxInit(getNumCols(), colsol);
01279   hasBeenOptimized_ = false;
01280 }
01281 
01292 void
01293 CouenneMINLPInterface::setRowPrice(const double * rowprice)
01294 {
01295   problem_->setDualsInit(getNumCols()*2 + getNumRows(), rowprice);
01296   hasBeenOptimized_ = false;
01297 }
01298 
01305 CoinWarmStart *
01306 CouenneMINLPInterface::getEmptyWarmStart () const
01307 {return app_->getEmptyWarmStart();}
01308 
01310 CoinWarmStart* 
01311 CouenneMINLPInterface::getWarmStart() const
01312 {
01313   if (exposeWarmStart_) {
01314     return internal_getWarmStart();;
01315   }
01316   else {
01317     return getEmptyWarmStart();
01318   }
01319 }
01322 bool 
01323 CouenneMINLPInterface::setWarmStart(const CoinWarmStart* ws)
01324 {
01325   if (exposeWarmStart_) {
01326     return internal_setWarmStart(ws);
01327   }
01328   else {
01329     return true;
01330   }
01331 }
01333 CoinWarmStart* 
01334 CouenneMINLPInterface::internal_getWarmStart() const
01335 {
01336   if (exposeWarmStart_ && warmstart_) {
01337     return warmstart_->clone();
01338   }
01339   else {
01340     return getEmptyWarmStart();
01341   }
01342 }
01345 bool 
01346 CouenneMINLPInterface::internal_setWarmStart(const CoinWarmStart* ws)
01347 {
01348   delete warmstart_;
01349   warmstart_ = NULL;
01350   hasBeenOptimized_ = false;
01351   if (exposeWarmStart_) {
01352     if (ws == NULL) {
01353       return true;
01354     }
01355     if(app_->warmStartIsValid(ws)) {
01356       warmstart_ = ws->clone();
01357       return true;
01358     }
01359     // See if it is anything else than the CoinWarmStartBasis that all others
01360     // derive from
01361     const CoinWarmStartPrimalDual* pdws =
01362       dynamic_cast<const CoinWarmStartPrimalDual*>(ws);
01363     if (pdws) {
01364       // Must be an IpoptWarmStart, since the others do not derive from this.
01365       // Accept it.
01366       warmstart_ = new IpoptWarmStart(*pdws);
01367       return true;
01368     }
01369     return false;
01370   }
01371   else {
01372     return true;
01373   }
01374 }
01375 
01377 void
01378 CouenneMINLPInterface::setContinuous(int index)
01379 {
01380   problem_->SetVariableType(index, TMINLP::CONTINUOUS);
01381   hasBeenOptimized_ = false;
01382 }
01384 void
01385 CouenneMINLPInterface::setInteger(int index)
01386 {
01387   problem_->SetVariableType(index, TMINLP::INTEGER);
01388   hasBeenOptimized_ = false;
01389 }
01390 
01392 double
01393 CouenneMINLPInterface::getObjValue() const
01394 {
01395   return problem_->obj_value();
01396 }
01397 
01398 //#############################################################################
01399 // Parameter related methods
01400 //#############################################################################
01401 
01402 bool
01403 CouenneMINLPInterface::setIntParam(OsiIntParam key, int value)
01404 {
01405   //  debugMessage("OsiCpxSolverInterface::setIntParam(%d, %d)\n", key, value);
01406 
01407   bool retval = false;
01408   switch (key) {
01409   case OsiMaxNumIteration:
01410     retval = false;
01411     break;
01412   case OsiMaxNumIterationHotStart:
01413     if( value >= 0 ) {
01414       retval = false;
01415     }
01416     else
01417       retval = false;
01418     break;
01419   case OsiLastIntParam:
01420     retval = false;
01421     break;
01422   default:
01423     retval = false;
01424     (*handler_)<< "Unhandled case in setIntParam\n"<<CoinMessageEol;
01425     break;
01426   }
01427   return retval;
01428 }
01429 
01430 //-----------------------------------------------------------------------------
01431 
01432 bool
01433 CouenneMINLPInterface::setDblParam(OsiDblParam key, double value)
01434 {
01435   //  debugMessage("CouenneMINLPInterface::setDblParam(%d, %g)\n", key, value);
01436 
01437   bool retval = false;
01438   switch (key) {
01439   case OsiDualObjectiveLimit:
01440     OsiDualObjectiveLimit_ = value;
01441     retval = true;
01442     break;
01443   case OsiPrimalObjectiveLimit:
01444     (*handler_)<<"Can not set primal objective limit parameter"<<CoinMessageEol;
01445     retval = false;
01446     break;
01447   case OsiDualTolerance:
01448     (*handler_)<<"Can not set dual tolerance parameter"<<CoinMessageEol;
01449     retval = false;
01450     break;
01451   case OsiPrimalTolerance:
01452     (*handler_)<<"Can not set primal tolerance parameter"<<CoinMessageEol;
01453     retval = false;
01454   case OsiObjOffset:
01455     retval = OsiSolverInterface::setDblParam(key,value);
01456     break;
01457   case OsiLastDblParam:
01458     retval = false;
01459     break;
01460   default:
01461     retval = false;
01462     (*handler_) << "Unhandled case in setDblParam"<<CoinMessageEol;
01463     break;
01464   }
01465   return retval;
01466 }
01467 
01468 
01469 //-----------------------------------------------------------------------------
01470 
01471 bool
01472 CouenneMINLPInterface::setStrParam(OsiStrParam key, const std::string & value)
01473 {
01474   //  debugMessage("CouenneMINLPInterface::setStrParam(%d, %s)\n", key, value.c_str());
01475 
01476   bool retval=false;
01477   switch (key) {
01478   case OsiProbName:
01479     OsiSolverInterface::setStrParam(key,value);
01480     return retval = true;
01481   case OsiSolverName:
01482     return false;
01483   case OsiLastStrParam:
01484     return false;
01485   }
01486   return false;
01487 }
01488 
01489 //-----------------------------------------------------------------------------
01490 
01491 bool
01492 CouenneMINLPInterface::getIntParam(OsiIntParam key, int& value) const
01493 {
01494   //  debugMessage("CouenneMINLPInterface::getIntParam(%d)\n", key);
01495 
01496   value = -COIN_INT_MAX; // Give a dummy value
01497   bool retval = false;
01498   switch (key) {
01499   case OsiMaxNumIteration:
01500     retval = false;
01501     break;
01502   case OsiMaxNumIterationHotStart:
01503     retval = false;
01504     break;
01505   case OsiLastIntParam:
01506     retval = false;
01507     break;
01508   default:
01509     retval = false;
01510     (*handler_) << "Unhandled case in setIntParam"<<CoinMessageEol;
01511   }
01512   return retval;
01513 }
01514 
01515 //-----------------------------------------------------------------------------
01516 
01517 bool
01518 CouenneMINLPInterface::getDblParam(OsiDblParam key, double& value) const
01519 {
01520   //  debugMessage("CouenneMINLPInterface::getDblParam(%d)\n", key);
01521 
01522   bool retval = false;
01523   switch (key) {
01524   case OsiDualObjectiveLimit:
01525     value = OsiDualObjectiveLimit_;
01526     retval = true;
01527     break;
01528   case OsiPrimalObjectiveLimit:
01529     value = getInfinity();
01530     retval = true;
01531     break;
01532   case OsiDualTolerance:
01533     retval = false;
01534     break;
01535   case OsiPrimalTolerance:
01536     options()->GetNumericValue("tol", value,"");
01537     //value = 1e-07;
01538     retval = true;
01539     break;
01540   case OsiObjOffset:
01541     retval = OsiSolverInterface::getDblParam(key, value);
01542     break;
01543   case OsiLastDblParam:
01544     retval = false;
01545     break;
01546   }
01547   return retval;
01548 }
01549 
01550 
01551 //-----------------------------------------------------------------------------
01552 
01553 bool
01554 CouenneMINLPInterface::getStrParam(OsiStrParam key, std::string & value) const
01555 {
01556   //  debugMessage("CouenneMINLPInterface::getStrParam(%d)\n", key);
01557 
01558   switch (key) {
01559   case OsiProbName:
01560     OsiSolverInterface::getStrParam(key, value);
01561     break;
01562   case OsiSolverName:
01563     value = "Ipopt";
01564     break;
01565   case OsiLastStrParam:
01566     return false;
01567   }
01568 
01569   return true;
01570 }
01571 
01572 void
01573 CouenneMINLPInterface::randomStartingPoint()
01574 {
01575   int numcols = getNumCols();
01576   const double * colLower = getColLower();
01577   const double * colUpper = getColUpper();
01578   double * sol = new double[numcols];
01579   const Number * x_init = problem_->x_init_user();
01580   const double* perturb_radius = NULL;
01581   if (randomGenerationType_ == perturb_suffix) {
01582     const TMINLP::PerturbInfo* pertubinfo = tminlp_->perturbInfo();
01583     if (pertubinfo) {
01584       perturb_radius = pertubinfo->GetPerturbationArray();
01585     }
01586     if (!perturb_radius) {
01587       throw SimpleError("Can't use perturb_radius if no radii are given.",
01588                         "randomStartingPoint");
01589     }
01590   }
01591   for(int i = 0 ; i < numcols ; i++) {
01592     int randomGenerationType = randomGenerationType_;
01593     if(x_init[i] < colLower[i] || x_init[i] > colUpper[i])
01594       randomGenerationType = uniform;
01595     if(randomGenerationType_ == uniform){
01596       double lower = std::min(-maxRandomRadius_,colUpper[i] - maxRandomRadius_);
01597       lower = std::max(colLower[i], lower);
01598       double upper = std::max(maxRandomRadius_,colLower[i] + maxRandomRadius_);
01599       upper = std::min(colUpper[i],upper);
01600       lower = std::min(upper,lower);
01601       upper = std::max(upper, lower);
01602       double interval = upper - lower;
01603       sol[i] = CoinDrand48()*(interval) + lower;}
01604     else if (randomGenerationType_ == perturb){
01605       const double lower = std::max(x_init[i] - max_perturbation_, colLower[i]);
01606       const double upper = std::min(x_init[i] + max_perturbation_, colUpper[i]);
01607       const double interval = upper - lower;
01608       sol[i]  = lower + CoinDrand48()*(interval);
01609     }
01610     else if (randomGenerationType_ == perturb_suffix){
01611       const double radius = perturb_radius[i];
01612       const double lower = std::max(x_init[i] - radius*max_perturbation_, colLower[i]);
01613       const double upper = std::min(x_init[i] + radius*max_perturbation_, colUpper[i]);
01614       const double interval = upper - lower;
01615       sol[i]  = lower + CoinDrand48()*(interval);
01616     }
01617   }
01618   app_->disableWarmStart();
01619   setColSolution(sol);
01620   delete [] sol;
01621 }
01622 
01623 
01624 
01626 int CouenneMINLPInterface::initializeJacobianArrays()
01627 {
01628   Index n, m, nnz_h_lag;
01629   TNLP::IndexStyleEnum index_style;
01630   tminlp_->get_nlp_info( n, m, nnz_jac, nnz_h_lag, index_style);
01631 
01632   if(jRow_ != NULL) delete jRow_;
01633   if(jCol_ != NULL) delete jCol_;
01634   if(jValues_ != NULL) delete jValues_;
01635 
01636   jRow_ = new Index[nnz_jac];
01637   jCol_ = new Index[nnz_jac];
01638   jValues_ = new Number[nnz_jac];
01639   tminlp_->eval_jac_g(n, NULL, 0, m, nnz_jac, jRow_, jCol_, NULL);
01640   if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
01641     {
01642       for(int i = 0 ; i < nnz_jac ; i++){
01643         jRow_[i]--;
01644         jCol_[i]--;
01645       }
01646     }
01647 
01648   if(constTypes_ != NULL) delete [] constTypes_;
01649   //  if(constTypesNum_ != NULL) delete [] constTypesNum_;
01650 
01651   constTypes_ = new TNLP::LinearityType[getNumRows()];
01652   tminlp_->get_constraints_linearity(getNumRows(), constTypes_);
01653   //  constTypesNum_ = new int[getNumRows()];
01654   for(int i = 0; i < getNumRows() ; i++) {
01655     if(constTypes_[i]==TNLP::NON_LINEAR) {
01656       //constTypesNum_[i] = 
01657       nNonLinear_++;
01658     }
01659   }
01660   return nnz_jac;
01661 }
01662 
01663 
01664 double 
01665 CouenneMINLPInterface::getConstraintsViolation(const double *x, double &obj)
01666 {
01667   int numcols = getNumCols();
01668   int numrows = getNumRows();
01669   double * g = new double[numrows];
01670   tminlp_->eval_g(numcols, x, 1, numrows, g);
01671   const double * rowLower = getRowLower();
01672   const double * rowUpper = getRowUpper();
01673 
01674 
01675   double norm = 0;
01676   for(int i = 0; i< numrows ; i++) {
01677     if(!constTypes_ || constTypes_[i] == TNLP::NON_LINEAR) {
01678       double rowViolation = 0;
01679       if(rowLower[i] > -1e10)
01680         rowViolation = std::max(0.,rowLower[i] - g[i]);
01681 
01682       if(rowUpper[i] < 1e10);
01683       rowViolation = std::max(rowViolation, g[i] - rowUpper[i]);
01684 
01685       norm = rowViolation > norm ? rowViolation : norm;
01686     }
01687   }
01688   tminlp_->eval_f(numcols, x, 1, obj);
01689   delete [] g;
01690   return norm;
01691 }
01692 
01694 double
01695 CouenneMINLPInterface::getNonLinearitiesViolation(const double *x, const double obj)
01696 {
01697   double f;
01698   double norm = getConstraintsViolation(x, f);
01699   assert((f - obj) > -1e-08);
01700   norm =  (f - obj) > norm ? f - obj : norm;
01701   return norm;
01702 }
01703 
01704 
01705 
01706 //A procedure to try to remove small coefficients in OA cuts (or make it non small
01707 static inline
01708 bool cleanNnz(double &value, double colLower, double colUpper,
01709               double rowLower, double rowUpper, double colsol,
01710               double & lb, double &ub, double tiny, double veryTiny)
01711 {
01712   if(fabs(value)>= tiny) return 1;
01713 
01714   if(fabs(value)<veryTiny) return 0;//Take the risk?
01715 
01716   //try and remove
01717   double infty = 1e20;
01718   bool colUpBounded = colUpper < 10000;
01719   bool colLoBounded = colLower > -10000;
01720   bool rowNotLoBounded =  rowLower <= - infty;
01721   bool rowNotUpBounded = rowUpper >= infty;
01722   bool pos =  value > 0;
01723 
01724   if(colLoBounded && pos && rowNotUpBounded) {
01725     lb += value * (colsol - colLower);
01726     return 0;
01727   }
01728   else
01729     if(colLoBounded && !pos && rowNotLoBounded) {
01730       ub += value * (colsol - colLower);
01731       return 0;
01732     }
01733     else
01734       if(colUpBounded && !pos && rowNotUpBounded) {
01735         lb += value * (colsol - colUpper);
01736         return 0;
01737       }
01738       else
01739         if(colUpBounded && pos && rowNotLoBounded) {
01740           ub += value * (colsol - colUpper);
01741           return 0;
01742         }
01743   //can not remove coefficient increase it to smallest non zero
01744   if(pos) value = tiny;
01745   else
01746     value = - tiny;
01747   return 1;
01748 }
01749 
01752 void
01753 CouenneMINLPInterface::getOuterApproximation(OsiCuts &cs, const double * x, bool getObj, const double * x2, double theta, bool global)
01754 {
01755 
01756   //printf ("this procedure has been called\n");
01757 
01758   int n,m, nnz_jac_g, nnz_h_lag;
01759   TNLP::IndexStyleEnum index_style;
01760   problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
01761   if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
01762     initializeJacobianArrays();
01763   assert(jRow_ != NULL);
01764   assert(jCol_ != NULL);
01765   double * g = new double[m];
01766 
01767   problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
01768   problem_to_optimize_->eval_g(n,x,1,m,g);
01769 
01770   //double *xx = new double [n];
01771 
01772   // x ---> xx 
01773   // 
01774   // 
01775   // 
01776 
01777 
01778 
01779 
01780   problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
01781   problem_to_optimize_->eval_g(n,x,1,m,g);
01782   //As jacobian is stored by cols fill OsiCuts with cuts
01783   CoinPackedVector * cuts = new CoinPackedVector[nNonLinear_ + 1];
01784   double * lb = new double[nNonLinear_ + 1];
01785   double * ub = new double[nNonLinear_ + 1];
01786 
01787   int * row2cutIdx = new int[m];//store correspondance between index of row and index of cut (some cuts are not generated for rows because linear, or not binding). -1 if constraint does not generate a cut, otherwise index in cuts.
01788   int numCuts = 0;
01789 
01790   const double * rowLower = getRowLower();
01791   const double * rowUpper = getRowUpper();
01792   const double * colLower = getColLower();
01793   const double * colUpper = getColUpper();
01794   const double * duals = getRowPrice() + 2 * n;
01795 
01796   for (int i=0; i< n; i++)
01797     printf ("variable %d: [%g,%g]\n", i, rowLower [i], rowUpper [i]);
01798 
01799   double infty = getInfinity();
01800   double nlp_infty = infty_;
01801   
01802   for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
01803     if(constTypes_[rowIdx] == TNLP::NON_LINEAR) {
01804 #if 0
01805       if(fabs(duals[rowIdx]) == 0.)
01806         {
01807           row2cutIdx[rowIdx] = -1;
01808 #ifdef NDEBUG
01809           (*handler_)<<"non binding constraint"<<CoinMessageEol;
01810 #endif
01811           continue;
01812         }
01813 #endif
01814       row2cutIdx[rowIdx] = numCuts;
01815       if(rowLower[rowIdx] > - nlp_infty)
01816         lb[numCuts] = rowLower[rowIdx] - g[rowIdx];
01817       else
01818         lb[numCuts] = - infty;
01819       if(rowUpper[rowIdx] < nlp_infty)
01820         ub[numCuts] = rowUpper[rowIdx] - g[rowIdx];
01821       else
01822         ub[numCuts] = infty;
01823       if(rowLower[rowIdx] > -infty && rowUpper[rowIdx] < infty)
01824         {
01825           if(duals[rowIdx] >= 0)// <= inequality
01826             lb[numCuts] = - infty;
01827           if(duals[rowIdx] <= 0)// >= inequality
01828             ub[numCuts] = infty;
01829         }
01830       
01831       numCuts++;
01832     }
01833     else
01834       row2cutIdx[rowIdx] = -1;
01835   }
01836 
01837 
01838   for(int i = 0 ; i < nnz_jac_g ; i++) {
01839     const int &rowIdx = jRow_[i];
01840     const int & cutIdx = row2cutIdx[ rowIdx ];
01841     if(cutIdx != -1) {
01842       const int &colIdx = jCol_[i];
01843       //"clean" coefficient
01844       if(cleanNnz(jValues_[i],colLower[colIdx], colUpper[colIdx],
01845                   rowLower[rowIdx], rowUpper[rowIdx],
01846                   x[colIdx],
01847                   lb[cutIdx],
01848                   ub[cutIdx], tiny_, veryTiny_)) {
01849         cuts[cutIdx].insert(colIdx,jValues_[i]);
01850         if(lb[cutIdx] > - infty)
01851           lb[cutIdx] += jValues_[i] * x[colIdx];
01852         if(ub[cutIdx] < infty)
01853           ub[cutIdx] += jValues_[i] * x[colIdx];
01854       }
01855     }
01856   }
01857 
01858   int * cut2rowIdx = NULL;
01859   if (IsValid(cutStrengthener_) || oaHandler_->logLevel() > 0) {
01860     cut2rowIdx = new int [numCuts];// Store correspondance between indices of cut and indices of rows. For each cut
01861     for(int rowIdx = 0 ; rowIdx < m ; rowIdx++){
01862       if(row2cutIdx[rowIdx] >= 0){
01863         cut2rowIdx[row2cutIdx[rowIdx]] = rowIdx;
01864       }
01865     }
01866   }
01867 
01868   for(int cutIdx = 0; cutIdx < numCuts; cutIdx++) {
01869     //Compute cut violation
01870     if(x2 != NULL) {
01871       double rhs = cuts[cutIdx].dotProduct(x2);
01872       double violation = 0.;
01873       violation = std::max(violation, rhs - ub[cutIdx]);
01874       violation = std::max(violation, lb[cutIdx] - rhs);
01875       if(violation < theta) {
01876         if(oaHandler_->logLevel() > 0)
01877           oaHandler_->message(CUT_NOT_VIOLATED_ENOUGH, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
01878         continue;}
01879       if(oaHandler_->logLevel() > 0)
01880         oaHandler_->message(VIOLATED_OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
01881     }
01882     else if (oaHandler_->logLevel() > 0)
01883       oaHandler_->message(OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<CoinMessageEol;
01884     OsiRowCut newCut;
01885     //    if(lb[i]>-1e20) assert (ub[i]>1e20);
01886 
01887     if (IsValid(cutStrengthener_)) {
01888       const int& rowIdx = cut2rowIdx[cutIdx];
01889       bool retval =
01890         cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
01891                                       GetRawPtr(problem_), rowIdx,
01892                                       cuts[cutIdx], lb[cutIdx], ub[cutIdx], g[rowIdx],
01893                                       rowLower[rowIdx], rowUpper[rowIdx],
01894                                       n, x, infty);
01895       if (!retval) {
01896         (*messageHandler()) << "error in cutStrengthener_->ComputeCuts\n";
01897         //exit(-2);
01898       }
01899     }
01900     if(global) {
01901       newCut.setGloballyValidAsInteger(1);
01902     }
01903     newCut.setEffectiveness(99.99e99);
01904     newCut.setLb(lb[cutIdx]);
01905     newCut.setUb(ub[cutIdx]);
01906     newCut.setRow(cuts[cutIdx]);
01907     //    CftValidator validator;
01908     //    validator(newCut);
01909     if(oaHandler_->logLevel()>2){
01910       oaHandler_->print(newCut);}
01911     cs.insert(newCut);
01912   }
01913 
01914   delete[] g;
01915   delete [] cuts;
01916   delete [] row2cutIdx;
01917   delete [] cut2rowIdx;
01918 
01919   if(getObj && ! problem_->hasLinearObjective()) { // Get the objective cuts
01920     double * obj = new double [n];
01921     problem_to_optimize_->eval_grad_f(n, x, 1,obj);
01922     double f;
01923     problem_to_optimize_->eval_f(n, x, 1, f);
01924 
01925     CoinPackedVector v;
01926     v.reserve(n);
01927     lb[nNonLinear_] = -f;
01928     ub[nNonLinear_] = -f;
01929     //double minCoeff = 1e50;
01930     for(int i = 0; i<n ; i++) {
01931       if(cleanNnz(obj[i],colLower[i], colUpper[i],
01932                   -getInfinity(), 0,
01933                   x[i],
01934                   lb[nNonLinear_],
01935                   ub[nNonLinear_],tiny_, 1e-15)) {
01936         //            minCoeff = std::min(fabs(obj[i]), minCoeff);
01937         v.insert(i,obj[i]);
01938         lb[nNonLinear_] += obj[i] * x[i];
01939         ub[nNonLinear_] += obj[i] * x[i];
01940       }
01941     }
01942     v.insert(n,-1);
01943     //Compute cut violation
01944     bool genCut = true;
01945     if(x2 != NULL) {
01946       double rhs = v.dotProduct(x2);
01947       double violation = std::max(0., rhs - ub[nNonLinear_]);
01948       if(violation < theta) genCut = false;
01949     }
01950     if(genCut) {
01951       if (IsValid(cutStrengthener_)) {
01952         lb[nNonLinear_] = -infty;
01953         bool retval =
01954           cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
01955                                         GetRawPtr(problem_), -1,
01956                                         v, lb[nNonLinear_], ub[nNonLinear_],
01957                                         ub[nNonLinear_], -infty, 0.,
01958                                         n, x, infty);
01959         if (!retval) {
01960           (*handler_)<< "error in cutStrengthener_->ComputeCuts"<<CoinMessageEol;
01961           //exit(-2);
01962         }
01963       }
01964       OsiRowCut newCut;
01965       if(global)
01966         newCut.setGloballyValidAsInteger(1);
01967       newCut.setEffectiveness(99.99e99);
01968       newCut.setRow(v);
01969       newCut.setLb(-COIN_DBL_MAX/*Infinity*/);
01970       newCut.setUb(ub[nNonLinear_]);
01971       //     CftValidator validator;
01972       //     validator(newCut);
01973       cs.insert(newCut);
01974     }
01975     delete [] obj;
01976   }
01977 
01978   delete []lb;
01979   delete[]ub;
01980 }
01981 
01982 
01983 
01986 void
01987 CouenneMINLPInterface::getConstraintOuterApproximation(OsiCuts &cs, int rowIdx, 
01988                                                        const double * x, 
01989                                                        const double * x2, bool global)
01990 {
01991   double g;
01992   int * indices = new int[getNumCols()];
01993   double * values = new double[getNumCols()];
01994   int nnz;
01995   problem_->eval_grad_gi(getNumCols(), x, 1, rowIdx, nnz, indices, values);
01996   problem_->eval_gi(getNumCols(),x,1, rowIdx, g);
01997 
01998   CoinPackedVector cut;
01999   double lb;
02000   double ub;
02001 
02002 
02003   const double rowLower = getRowLower()[rowIdx];
02004   const double rowUpper = getRowUpper()[rowIdx];
02005   const double * colLower = getColLower();
02006   const double * colUpper = getColUpper();
02007   const double dual = (getRowPrice() + 2 * getNumCols())[rowIdx];
02008   double infty = getInfinity();
02009   double nlp_infty = infty_;
02010   
02011   if(rowLower > - nlp_infty)
02012     lb = rowLower - g;
02013   else
02014     lb = - infty;
02015   if(rowUpper < nlp_infty)
02016     ub = rowUpper - g;
02017   else
02018     ub = infty;
02019   if(rowLower > -infty && rowUpper < infty)
02020     {
02021       if(dual >= 0)// <= inequality
02022         lb = - infty;
02023       if(dual <= 0)// >= inequality
02024         ub = infty;
02025     }
02026 
02027   for(int i = 0 ; i < nnz; i++) {
02028     const int &colIdx = indices[i];
02029     //"clean" coefficient
02030     if(cleanNnz(values[i],colLower[colIdx], colUpper[colIdx],
02031                 rowLower, rowUpper,
02032                 x[colIdx],
02033                 lb,
02034                 ub, tiny_, veryTiny_)) {
02035       cut.insert(colIdx,values[i]);
02036       if(lb > - infty)
02037         lb += values[i] * x[colIdx];
02038       if(ub < infty)
02039         ub += values[i] * x[colIdx];
02040     }
02041   }
02042 
02043   OsiRowCut newCut;
02044 
02045   if(global) {
02046     newCut.setGloballyValidAsInteger(1);
02047   }
02048   newCut.setEffectiveness(99.99e99);
02049   newCut.setLb(lb);
02050   newCut.setUb(ub);
02051   newCut.setRow(cut);
02052   cs.insert(newCut);
02053 
02054   delete [] indices;
02055   delete [] values;
02056 }
02057 
02058 void
02059 CouenneMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
02060                                                   double a, double s, int L){
02061   if(! IsValid(feasibilityProblem_)) {
02062     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02063   }
02064   feasibilityProblem_->set_use_feasibility_pump_objective(true);
02065   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02066   feasibilityProblem_->setLambda(a);
02067   feasibilityProblem_->setSigma(s);
02068   feasibilityProblem_->setNorm(L);
02069   feasibilityProblem_->set_use_cutoff_constraint(false);
02070   feasibilityProblem_->set_use_local_branching_constraint(false);  
02071   problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
02072   feasibility_mode_ = true;
02073 }
02074 
02075 void
02076 CouenneMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
02077                                                   double rhs_local_branching_constraint){
02078   if(! IsValid(feasibilityProblem_)) {
02079     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02080   }
02081   feasibilityProblem_->set_use_feasibility_pump_objective(false);
02082   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02083   feasibilityProblem_->set_use_cutoff_constraint(false);
02084   feasibilityProblem_->set_use_local_branching_constraint(true);  
02085   feasibilityProblem_->set_rhs_local_branching_constraint(rhs_local_branching_constraint);  
02086   problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
02087   feasibility_mode_ = true;
02088 }
02089 
02090 void
02091 CouenneMINLPInterface::switchToOriginalProblem(){
02092   problem_to_optimize_ = GetRawPtr(problem_);
02093   feasibility_mode_ = false;
02094 }
02095 
02096 double
02097 CouenneMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
02098                                                double a, double s, int L)
02099 {
02100   if(! IsValid(feasibilityProblem_)) {
02101     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02102   }
02103   feasibilityProblem_->set_use_feasibility_pump_objective(true);
02104   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02105   feasibilityProblem_->setLambda(a);
02106   feasibilityProblem_->setSigma(s);
02107   feasibilityProblem_->setNorm(L);
02108   feasibilityProblem_->set_use_cutoff_constraint(false);
02109   feasibilityProblem_->set_use_local_branching_constraint(false);  
02110   nCallOptimizeTNLP_++;
02111   totalNlpSolveTime_-=CoinCpuTime();
02112   SmartPtr<TNLPSolver> app2 = app_->clone();
02113   app2->options()->SetIntegerValue("print_level", (Index) 0);
02114   optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
02115   totalNlpSolveTime_+=CoinCpuTime();
02116   hasBeenOptimized_=true;
02117   return getObjValue();
02118 }
02119 
02120 double
02121 CouenneMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
02122                                                int L, double cutoff)
02123 {
02124   if(! IsValid(feasibilityProblem_)) {
02125     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02126   }
02127   feasibilityProblem_->set_use_feasibility_pump_objective(true);
02128   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02129   feasibilityProblem_->setLambda(1.0);
02130   feasibilityProblem_->setSigma(0.0);
02131   feasibilityProblem_->setNorm(L);
02132   feasibilityProblem_->set_use_cutoff_constraint(true);
02133   feasibilityProblem_->set_cutoff(cutoff);
02134   feasibilityProblem_->set_use_local_branching_constraint(false);  
02135   nCallOptimizeTNLP_++;
02136   totalNlpSolveTime_-=CoinCpuTime();
02137   SmartPtr<TNLPSolver> app2 = app_->clone();
02138   app2->options()->SetIntegerValue("print_level", (Index) 0);
02139   optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
02140   totalNlpSolveTime_+=CoinCpuTime();
02141   hasBeenOptimized_=true;
02142   return getObjValue();
02143 }
02144 
02145 double
02146 CouenneMINLPInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs, bool addOnlyViolated, bool global)
02147 {
02148   double ret_val = solveFeasibilityProblem(n, x_bar, inds, 1, 0, 2);
02149   getOuterApproximation(cs, getColSolution(), 0, (addOnlyViolated? x_bar:NULL)
02150                         , global);
02151   return ret_val;
02152 }
02153 
02154 
02155 static bool WarnedForNonConvexOa=false;
02156 
02157 void
02158 CouenneMINLPInterface::extractLinearRelaxation(OsiSolverInterface &si, 
02159                                                const double * x, bool getObj)
02160 {
02161   double * rowLow = NULL;
02162   double * rowUp = NULL;
02163 
02164   int n;
02165   int m;
02166   int nnz_jac_g;
02167   int nnz_h_lag;
02168   TNLP::IndexStyleEnum index_style;
02169   //Get problem information
02170   problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
02171 
02172   //if not allocated allocate spaced for stroring jacobian
02173   if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
02174     initializeJacobianArrays();
02175 
02176   //get Jacobian
02177   problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
02178 
02179 
02180   double *g = new double[m];
02181   problem_to_optimize_->eval_g(n, x, 1, m, g);
02182 
02183   rowLow = new double[m];
02184   rowUp = new double[m];
02185   int * nonBindings = new int[m];//store non binding constraints (which are to be removed from OA)
02186   int numNonBindings = 0;
02187   const double * rowLower = getRowLower();
02188   const double * rowUpper = getRowUpper();
02189   const double * colLower = getColLower();
02190   const double * colUpper = getColUpper();
02191   const double * duals = getRowPrice() + 2*n;
02192   assert(m==getNumRows());
02193   double infty = si.getInfinity();
02194   double nlp_infty = infty_;
02195   
02196   for(int i = 0 ; i < m ; i++) {
02197     if(constTypes_[i] == TNLP::NON_LINEAR) {
02198       //If constraint is range not binding prepare to remove it
02199       if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty && fabs(duals[i]) == 0.)
02200         {
02201           nonBindings[numNonBindings++] = i;
02202           continue;
02203         }
02204       else
02205         if(rowLower[i] > - nlp_infty){
02206           rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
02207           if(! WarnedForNonConvexOa && rowUpper[i] < nlp_infty){
02208             messageHandler()->message(WARNING_NON_CONVEX_OA, messages_)<<CoinMessageEol;
02209             WarnedForNonConvexOa = true;
02210           }
02211         }
02212         else
02213           rowLow[i] = - infty;
02214       if(rowUpper[i] < nlp_infty)
02215         rowUp[i] =  (rowUpper[i] - g[i]) + 1e-07;
02216       else
02217         rowUp[i] = infty;
02218       
02219       //If equality or ranged constraint only add one side by looking at sign of dual multiplier
02220       if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty)
02221         {
02222           if(duals[i] >= 0.)// <= inequality
02223             rowLow[i] = - infty;
02224           if(duals[i] <= 0.)// >= inequality
02225             rowUp[i] = infty;
02226         }
02227     }
02228     else {
02229       if(rowLower[i] > -nlp_infty){
02230         //   printf("Lower %g ", rowLower[i]);
02231         rowLow[i] = (rowLower[i] - g[i]);
02232       }
02233       else
02234         rowLow[i] = - infty;
02235       if(rowUpper[i] < nlp_infty){
02236         //   printf("Upper %g ", rowUpper[i]);
02237         rowUp[i] =  (rowUpper[i] - g[i]);
02238       }
02239       else
02240         rowUp[i] = infty;
02241     }
02242   }
02243 
02244   
02245   
02246   //Then convert everything to a CoinPackedMatrix
02247   //Go through values, clean coefficients and fix bounds
02248   for(int i = 0 ; i < nnz_jac_g ; i++) {
02249     if(constTypes_[jRow_[i]] != TNLP::LINEAR){//For linear just copy is fine.
02250       if(//For other clean tinys
02251          cleanNnz(jValues_[i],colLower[jCol_[i]], colUpper[jCol_[i]],
02252                   rowLower[jRow_[i]], rowUpper[jRow_[i]],
02253                   x[jCol_[i]],
02254                   rowLow[jRow_[i]],
02255                   rowUp[jRow_[i]], tiny_, veryTiny_)) {      
02256         rowLow[jRow_[i]] += jValues_[i] * x[jCol_ [i]];
02257         rowUp[jRow_[i]] += jValues_[i] *x[jCol_[i]];
02258       }
02259     }
02260     else {
02261       double value = jValues_[i] * getColSolution()[jCol_[i]];
02262       rowLow[jRow_[i]] += value;
02263       rowUp[jRow_[i]] += value;
02264     } 
02265   }
02266   CoinPackedMatrix mat(true, jRow_, jCol_, jValues_, nnz_jac_g);
02267   mat.setDimensions(m,n); // In case matrix was empty, this should be enough
02268   
02269   //remove non-bindings equality constraints
02270   mat.deleteRows(numNonBindings, nonBindings);
02271 
02272   int numcols=getNumCols();
02273   double *obj = new double[numcols];
02274   for(int i = 0 ; i < numcols ; i++)
02275     obj[i] = 0.;
02276   
02277   
02278   si.loadProblem(mat, getColLower(), getColUpper(), obj, rowLow, rowUp);
02279   delete [] rowLow;
02280   delete [] rowUp;
02281   delete [] nonBindings;
02282   delete [] g;
02283   for(int i = 0 ; i < getNumCols() ; i++) {
02284     if(isInteger(i))
02285       si.setInteger(i);
02286   }
02287   if(getObj) {
02288     bool addObjVar = false;
02289     if(problem_->hasLinearObjective()){
02290       //Might be in trouble if objective has a constant part
02291       // for now just check that f(0) = 0. 
02292       // If it is not adding a constant term does not seem supported by Osi
02293       // for now
02294       double zero;
02295       problem_to_optimize_->eval_f(n, obj, 1, zero);
02296       si.setDblParam(OsiObjOffset, -zero);
02297       //if(fabs(zero - 0) > 1e-10)
02298       //addObjVar = true;
02299       //else { 
02300       //Copy the linear objective and don't create a dummy variable.
02301       problem_to_optimize_->eval_grad_f(n, x, 1,obj);
02302       si.setObjective(obj);
02303       //}
02304     }
02305     else {
02306       addObjVar = true;
02307     }
02308 
02309     if(addObjVar){
02310       //add variable alpha
02311       //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
02312       CoinPackedVector a;
02313       si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
02314   
02315       // Now get the objective cuts
02316       // get the gradient, pack it and add the cut
02317       problem_to_optimize_->eval_grad_f(n, x, 1,obj);
02318       double ub;
02319       problem_to_optimize_->eval_f(n, x, 1, ub);
02320       ub*=-1;
02321       double lb = -1e300;
02322       CoinPackedVector objCut;
02323       CoinPackedVector * v = &objCut;
02324       v->reserve(n);
02325       for(int i = 0; i<n ; i++) {
02326         if(nnz_jac_g)
02327           {
02328             if(cleanNnz(obj[i],colLower[i], colUpper[i],
02329                         -getInfinity(), 0,
02330                         x[i],
02331                         lb,
02332                         ub, tiny_, veryTiny_)) {
02333               v->insert(i,obj[i]);
02334               lb += obj[i] * x[i];
02335               ub += obj[i] * x[i];
02336             }
02337           }
02338         else //Unconstrained problem can not put clean coefficient
02339           {
02340             if(cleanNnz(obj[i],colLower[i], colUpper[i],
02341                         -getInfinity(), 0,
02342                         x[i],
02343                         lb,
02344                         ub, 1e-03, 1e-08)) {
02345               v->insert(i,obj[i]);
02346               lb += obj[i] * x[i];
02347               ub += obj[i] * x[i];
02348             }
02349           }
02350       }
02351       v->insert(n,-1);
02352       si.addRow(objCut, lb, ub);
02353     }
02354   }
02355   //  si.writeMpsNative("OA.mps",NULL, NULL, 1);
02356   delete [] obj;
02357 }
02358 
02360 void 
02361 CouenneMINLPInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
02362 { 
02363   if(numberCuts)
02364     freeCachedRowRim();
02365   const OsiRowCut ** cutsPtrs = new const OsiRowCut*[numberCuts];
02366   for(int i = 0 ; i < numberCuts ; i++)
02367     {
02368       cutsPtrs[i] = &cuts[i];
02369     }
02370   problem_->addCuts(numberCuts, cutsPtrs);
02371   delete [] cutsPtrs;
02372 }
02373 
02374 void
02375 CouenneMINLPInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
02376                                            const char * whereFrom)
02377 {
02378   totalNlpSolveTime_-=CoinCpuTime();
02379   if(warmStarted)
02380     optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
02381   else
02382     optimizationStatus_ = app_->OptimizeTNLP(GetRawPtr(problem_to_optimize_));
02383   totalNlpSolveTime_+=CoinCpuTime();
02384   nCallOptimizeTNLP_++;
02385   hasBeenOptimized_ = true;
02386  
02387  
02388   //Options should have been printed if not done already turn off Ipopt output
02389   if(!hasPrintedOptions) {
02390     hasPrintedOptions = 1;
02391     //app_->Options()->SetIntegerValue("print_level",0, true, true);
02392     app_->options()->SetStringValue("print_user_options","no", false, true);
02393   }
02394   
02395   bool otherDisagree = false ;
02396 #if 0
02397   if(optimizationStatus_ == TNLPSolver::notEnoughFreedom)//Too few degrees of freedom
02398     {
02399       (*messageHandler())<<"Too few degrees of freedom...."<<CoinMessageEol;
02400       int numrows = getNumRows();
02401       int numcols = getNumCols();
02402     
02403       const double * colLower = getColLower();
02404       const double * colUpper = getColUpper();
02405     
02406     
02407       const double * rowLower = getRowLower();
02408       const double * rowUpper = getRowUpper();
02409     
02410       int numberFixed = 0;
02411       for(int i = 0 ; i < numcols ; i++)
02412         {
02413           if(colUpper[i] - colLower[i] <= INT_BIAS)
02414             {
02415               numberFixed++;
02416             }
02417         }
02418       int numberEqualities = 0;
02419       for(int i = 0 ; i < numrows ; i++)
02420         {
02421           if(rowUpper[i] - rowLower[i] <= INT_BIAS)
02422             {
02423               numberEqualities++;
02424             }     
02425         }
02426       if(numcols - numberFixed > numberEqualities || numcols < numberEqualities)
02427         {
02428           std::string probName;
02429           getStrParam(OsiProbName, probName);
02430           throw newUnsolvedError(app_->errorCode(), problem_, probName);
02431         }
02432       double * saveColLow = CoinCopyOfArray(getColLower(), getNumCols());
02433       double * saveColUp = CoinCopyOfArray(getColUpper(), getNumCols());
02434     
02435       for(int i = 0; i < numcols && numcols - numberFixed <= numberEqualities ; i++)
02436         {
02437           if(colUpper[i] - colLower[i] <= INT_BIAS)
02438             {
02439               setColLower(i, saveColLow[i]-1e-06);
02440               setColUpper(i, saveColLow[i]+1e-06);
02441               numberFixed--;
02442             }
02443         }
02444       solveAndCheckErrors(warmStarted, throwOnFailure, whereFrom);
02445       //restore
02446       for(int i = 0; i < numcols && numcols - numberFixed < numrows ; i++)
02447         {
02448           problem_->SetVariableLowerBound(i,saveColLow[i]);
02449           problem_->SetVariableUpperBound(i,saveColUp[i]);
02450         }
02451       delete [] saveColLow;
02452       delete [] saveColUp;
02453       return;
02454     }
02455   else 
02456 #endif
02457     if(!app_->isRecoverable(optimizationStatus_))//Solver failed and the error can not be recovered, throw it
02458       {
02459         std::string probName;
02460         getStrParam(OsiProbName, probName);
02461         throw newUnsolvedError(app_->errorCode(), problem_, probName);
02462       }
02463     else if(testOthers_ && !app_->isError(optimizationStatus_)){
02464       Ipopt::SmartPtr<TMINLP2TNLP> problem_copy = problem_->clone();
02465       //Try other solvers and see if they agree
02466       int f =1;
02467       for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
02468           i != debug_apps_.end() ; i++){
02469         TNLPSolver::ReturnStatus otherStatus = (*i)->OptimizeTNLP(GetRawPtr(problem_copy));
02470         messageHandler()->message(LOG_LINE, messages_)
02471           <<'d'<<f++<<statusAsString(otherStatus)<<problem_copy->obj_value()
02472           <<(*i)->IterationCount()<<(*i)->CPUTime()<<CoinMessageEol;
02473         if(!(*i)->isError(otherStatus)){
02474           CoinRelFltEq eq(1e-05);
02475           if(otherStatus != optimizationStatus_){
02476             otherDisagree = true;
02477             messageHandler()->message(SOLVER_DISAGREE_STATUS, messages_)
02478               <<app_->solverName()<<statusAsString()
02479               <<(*i)->solverName()<<statusAsString(otherStatus)<<CoinMessageEol; 
02480           }
02481           else if(isProvenOptimal() && !eq(problem_->obj_value(),problem_copy->obj_value()))
02482             {
02483               otherDisagree = true;
02484               messageHandler()->message(SOLVER_DISAGREE_VALUE, messages_)
02485                 <<app_->solverName()<<problem_->obj_value()
02486                 <<(*i)->solverName()<<problem_copy->obj_value()<<CoinMessageEol; 
02487             }
02488         }
02489       }
02490     }
02491   try{
02492     totalIterations_ += app_->IterationCount();
02493   }
02494   catch(SimpleError &E)
02495     {
02496       if (throwOnFailure)//something failed throw
02497         {
02498           throw SimpleError("No statistics available from Ipopt",whereFrom);
02499         }
02500       else {
02501         return;
02502       }
02503     }
02504   if(problem_->hasUpperBoundingObjective()){//Check if solution is integer and recompute objective value using alternative objective function
02505     const double * sol = getColSolution();
02506     bool integerSol = true;
02507     double intTol = 1e-08;
02508     if(objects()){
02509       int nObjects = numberObjects();
02510       OsiObject ** object = objects();
02511       for(int i = 0 ; i< nObjects ; i++){
02512         int dummy;
02513         if(object[i]->infeasibility(this,dummy) > intTol)
02514           {
02515             integerSol=false;
02516             break;
02517           }
02518       }
02519     }
02520     else{//Only works for integer constraints
02521       int numcols = getNumCols();
02522       for(int i = 0 ; i < numcols ; i++){
02523         if(isInteger(i) || isBinary(i)){
02524           if(fabs(sol[i] - floor(sol[i]+0.5)) > intTol){
02525             integerSol = false;
02526             break;
02527           }
02528         }
02529       }
02530     }
02531     if(integerSol&&isProvenOptimal()){
02532       double help= problem_->evaluateUpperBoundingFunction(sol);
02533      
02534 
02535       OsiAuxInfo * auxInfo = getAuxiliaryInfo();
02536       Bonmin::AuxInfo * bonInfo = dynamic_cast<Bonmin::AuxInfo *>(auxInfo);
02537       if(bonInfo!=0)
02538         {
02539         
02540           if(help<bonInfo->bestObj2())
02541             {
02542               bonInfo->setBestObj2(help);
02543               bonInfo->setBestSolution2(getNumCols(), const_cast<double *>(getColSolution()));
02544 
02545               messageHandler()->message(ALTERNATE_OBJECTIVE, messages_)
02546                 <<help<<CoinMessageEol;
02547             }
02548         }
02549       else {
02550         printf("\nWARNING: the algorithm selected does not consider the second objective function\n");
02551       }
02552     }
02553   }
02554   messageHandler()->message(IPOPT_SUMMARY, messages_)
02555     <<whereFrom<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
02556   
02557   if((nCallOptimizeTNLP_ % 20) == 1)
02558     messageHandler()->message(LOG_HEAD, messages_)<<CoinMessageEol;
02559   
02560   
02561   if ( (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) ||
02562        ( otherDisagree )){
02563     messageHandler()->message(SUSPECT_PROBLEM,
02564                               messages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
02565     std::string subProbName;
02566     getStrParam(OsiProbName, subProbName);
02567     std::ostringstream os;
02568     os<<"_"<<nCallOptimizeTNLP_;
02569     subProbName+=os.str();
02570     problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
02571   }
02572   
02573 }
02574 
02576 // Solve Methods                                                  //
02579 void CouenneMINLPInterface::initialSolve()
02580 {
02581   assert(IsValid(app_));
02582   assert(IsValid(problem_));
02583 
02584   // Discard warmstart_ if we had one
02585   delete warmstart_;
02586   warmstart_ = NULL;
02587   
02588   if(!hasPrintedOptions) {
02589     int printOptions;
02590     app_->options()->GetEnumValue("print_user_options",printOptions,"bonmin.");
02591     if(printOptions)
02592       app_->options()->SetStringValue("print_user_options","yes",true,true);
02593   }
02594   if(exposeWarmStart_)
02595     app_->disableWarmStart(); 
02596   solveAndCheckErrors(0,1,"initialSolve");
02597   
02598   //Options should have been printed if not done already turn off Ipopt output
02599   if(!hasPrintedOptions) {
02600     hasPrintedOptions = 1;
02601     app_->options()->SetStringValue("print_user_options","no");
02602     app_->options()->SetIntegerValue("print_level",0);
02603   }
02604   
02605   messageHandler()->message(LOG_FIRST_LINE, messages_)<<nCallOptimizeTNLP_
02606                                                       <<statusAsString()
02607                                                       <<getObjValue()
02608                                                       <<app_->IterationCount()
02609                                                       <<app_->CPUTime()
02610                                                       <<CoinMessageEol;
02611   
02612   int numRetry = firstSolve_ ? numRetryInitial_ : numRetryResolve_;
02613   if(isAbandoned()) {
02614     resolveForRobustness(numRetryUnsolved_);
02615   }
02616   else if(numRetry)
02617     {
02618       resolveForCost(numRetry, numRetryInitial_ > 0);
02620       numRetryInitial_ = 0;
02621     }
02622   firstSolve_ = false;
02623 
02624   // if warmstart_ is not empty then had to use resolveFor... and that created
02625   // the warmstart at the end, and we have nothing to do here. Otherwise...
02626   if (! warmstart_ && ! isAbandoned()) {
02627     if (exposeWarmStart_) {
02628       warmstart_ = app_->getWarmStart(problem_);
02629     }
02630   }
02631 }
02632 
02635 void
02636 CouenneMINLPInterface::resolve()
02637 {
02638   assert(IsValid(app_));
02639   assert(IsValid(problem_));
02640   
02641   int has_warmstart = warmstart_ == NULL ? 0 : 1;
02642   if(warmstart_ == NULL) has_warmstart = 0;
02643   else if(!app_->warmStartIsValid(warmstart_)) has_warmstart = 1;
02644   else has_warmstart = 2;
02645   if (has_warmstart < 2) {
02646     // empty (or unrecognized) warmstart
02647     initialSolve();
02648     return;
02649   }
02650   app_->setWarmStart(warmstart_, problem_);
02651   delete warmstart_;
02652   warmstart_ = NULL;
02653 
02654   if (INT_BIAS > 0.) {
02655     app_->options()->SetStringValue("warm_start_same_structure", "yes");
02656   }
02657   else {
02658     app_->options()->SetStringValue("warm_start_same_structure", "no");
02659   }
02660 
02661   if(problem_->duals_init() != NULL)
02662     app_->enableWarmStart();
02663   else app_->disableWarmStart();
02664   solveAndCheckErrors(1,1,"resolve");
02665   
02666   messageHandler()->message(LOG_FIRST_LINE, messages_)<<nCallOptimizeTNLP_
02667                                                       <<statusAsString()
02668                                                       <<getObjValue()
02669                                                       <<app_->IterationCount()
02670                                                       <<app_->CPUTime()<<CoinMessageEol;
02671   
02672   if(isAbandoned()) {
02673     resolveForRobustness(numRetryUnsolved_);
02674   }
02675   else if(numRetryResolve_ ||
02676           (numRetryInfeasibles_ && isProvenPrimalInfeasible() ))
02677     resolveForCost(std::max(numRetryResolve_, numRetryInfeasibles_), 0);
02678 
02679   // if warmstart_ is not empty then had to use resolveFor... and that created
02680   // the warmstart at the end, and we have nothing to do here. Otherwise...
02681   if (! warmstart_ && ! isAbandoned()) {
02682     if (exposeWarmStart_) {
02683       warmstart_ = app_->getWarmStart(problem_);
02684     }
02685   }
02686 }
02687 
02688 
02690 // Methods returning info on how the solution process terminated  //
02693 bool CouenneMINLPInterface::isAbandoned() const
02694 {
02695   return (
02696           (optimizationStatus_==TNLPSolver::iterationLimit)||
02697           (optimizationStatus_==TNLPSolver::computationError)||
02698           (optimizationStatus_==TNLPSolver::illDefinedProblem)||
02699           (optimizationStatus_==TNLPSolver::illegalOption)||
02700           (optimizationStatus_==TNLPSolver::externalException)|
02701           (optimizationStatus_==TNLPSolver::exception)
02702           );
02703 }
02704 
02706 bool CouenneMINLPInterface::isProvenOptimal() const
02707 {
02708   return (optimizationStatus_==TNLPSolver::solvedOptimal) ||
02709     (optimizationStatus_==TNLPSolver::solvedOptimalTol);
02710 }
02712 bool CouenneMINLPInterface::isProvenPrimalInfeasible() const
02713 {
02714   return (optimizationStatus_ == TNLPSolver::provenInfeasible);
02715 }
02717 bool CouenneMINLPInterface::isProvenDualInfeasible() const
02718 {
02719   return (optimizationStatus_ == TNLPSolver::unbounded);
02720 }
02722 bool CouenneMINLPInterface::isPrimalObjectiveLimitReached() const
02723 {
02724   (*handler_)<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<CoinMessageEol;
02725   return 0;
02726 }
02728 bool CouenneMINLPInterface::isDualObjectiveLimitReached() const
02729 {
02730   //  (*messageHandler_)<<"Warning : isDualObjectiveLimitReached not implemented yet"<<CoinMessageEol;
02731   return (optimizationStatus_==TNLPSolver::unbounded);
02732 
02733 }
02735 bool CouenneMINLPInterface::isIterationLimitReached() const
02736 {
02737   return (optimizationStatus_==TNLPSolver::iterationLimit);
02738 }
02739 
02740 void
02741 CouenneMINLPInterface::extractInterfaceParams()
02742 {
02743   if (IsValid(app_)) {
02744     int logLevel;
02745     app_->options()->GetIntegerValue("nlp_log_level", logLevel,"bonmin.");
02746     messageHandler()->setLogLevel(logLevel);
02747 
02748 #ifdef COIN_HAS_FILTERSQP
02749     FilterSolver * filter = dynamic_cast<FilterSolver *>(GetRawPtr(app_));
02750 
02751     bool is_given =
02752 #endif
02753       app_->options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,"bonmin.");
02754 
02755 #ifdef COIN_HAS_FILTERSQP
02756     if(filter && !is_given){
02757       // overwriting default for filter
02758       maxRandomRadius_ = 10.;
02759     }
02760 #endif
02761    
02762     int oaCgLogLevel = 0;
02763     app_->options()->GetIntegerValue("oa_cuts_log_level", oaCgLogLevel,"bonmin.");
02764     oaHandler_->setLogLevel(oaCgLogLevel); 
02765     
02766     int exposeWs = false;
02767     app_->options()->GetEnumValue("warm_start", exposeWs, "bonmin.");
02768     setExposeWarmStart(exposeWs > 0);
02769     
02770     app_->options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,"bonmin.");
02771     app_->options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,"bonmin.");
02772     app_->options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,"bonmin.");
02773     app_->options()->GetIntegerValue("num_resolve_at_infeasibles", numRetryInfeasibles_,"bonmin.");
02774     app_->options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,"bonmin.");
02775     app_->options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,"bonmin.");
02776     app_->options()->GetNumericValue
02777       ("warm_start_bound_frac" ,pushValue_,"bonmin.");
02778     app_->options()->GetNumericValue("tiny_element",tiny_,"bonmin.");
02779     app_->options()->GetNumericValue("very_tiny_element",veryTiny_,"bonmin.");
02780     app_->options()->GetNumericValue("random_point_perturbation_interval",max_perturbation_,"bonmin.");
02781     app_->options()->GetEnumValue("random_point_type",randomGenerationType_,"bonmin.");
02782     int cut_strengthening_type;
02783     app_->options()->GetEnumValue("cut_strengthening_type", cut_strengthening_type,"bonmin.");
02784 
02785     if (cut_strengthening_type != CS_None) {
02786       // TNLP solver to be used in the cut strengthener
02787       cutStrengthener_ = new CutStrengthener(app_->clone(), app_->options());
02788     }
02789   }
02790 }
02791 
02792 void
02793 CouenneMINLPInterface::SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)
02794 {
02795   strong_branching_solver_ = strong_branching_solver;
02796 }
02797 
02798 //#define STRONG_COMPARE
02799 #ifdef STRONG_COMPARE
02800 static double objorig;
02801 #endif
02802 
02803 void
02804 CouenneMINLPInterface::markHotStart()
02805 {
02806   if (IsValid(strong_branching_solver_)) {
02807 #ifdef STRONG_COMPARE
02808     // AWDEBUG
02809     OsiSolverInterface::markHotStart();
02810     objorig = getObjValue();
02811 #endif
02812     optimizationStatusBeforeHotStart_ = optimizationStatus_;
02813     strong_branching_solver_->markHotStart(this);
02814   }
02815   else {
02816     // Default Implementation
02817     OsiSolverInterface::markHotStart();
02818   }
02819 }
02820 
02821 void
02822 CouenneMINLPInterface::solveFromHotStart()
02823 {
02824   if (IsValid(strong_branching_solver_)) {
02825 #ifdef STRONG_COMPARE
02826     // AWDEBUG
02827     OsiSolverInterface::solveFromHotStart();
02828     double obj_nlp = getObjValue() - objorig;
02829 #endif
02830     optimizationStatus_ = strong_branching_solver_->solveFromHotStart(this);
02831     hasBeenOptimized_ = true;
02832 #ifdef STRONG_COMPARE
02833     double obj_other = getObjValue() - objorig;
02834     printf("AWDEBUG: Strong Branching results: NLP = %15.8e Other = %15.8e\n",
02835            obj_nlp, obj_other);
02836 #endif
02837   }
02838   else {
02839     // Default Implementation
02840     OsiSolverInterface::solveFromHotStart();
02841   }
02842 }
02843 
02844 void
02845 CouenneMINLPInterface::unmarkHotStart()
02846 {
02847   if (IsValid(strong_branching_solver_)) {
02848 #ifdef STRONG_COMPARE
02849     // AWDEBUG
02850     OsiSolverInterface::unmarkHotStart();
02851 #endif
02852     strong_branching_solver_->unmarkHotStart(this);
02853     optimizationStatus_ = optimizationStatusBeforeHotStart_;
02854   }
02855   else {
02856     // Default Implementation
02857     OsiSolverInterface::unmarkHotStart();
02858   }
02859 }
02860 
02861 const double * CouenneMINLPInterface::getObjCoefficients() const
02862 {
02863   const int n = getNumCols();
02864   delete [] obj_;
02865   obj_ = NULL;
02866   obj_ = new double[n];
02867 
02868   bool new_x = true;
02869   const double* x_sol = problem_->x_sol();
02870   bool retval = problem_->eval_grad_f(n, x_sol, new_x, obj_);
02871   
02872   if (!retval) {
02873     // Let's see if that happens - it will cause a crash
02874     fprintf(stderr, "ERROR WHILE EVALUATING GRAD_F in CouenneMINLPInterface::getObjCoefficients()\n");
02875     delete [] obj_;
02876     obj_ = NULL;
02877   }
02878 
02879   return obj_;
02880 }
02881 
02882 
02883 }
02885 #endif

Generated on Wed Nov 30 03:04:05 2011 by  doxygen 1.4.7