/home/coin/SVN-release/OS-1.1.1/Bonmin/src/Interfaces/BonOsiTMINLPInterface.cpp

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

Generated on Tue Sep 30 03:01:23 2008 by  doxygen 1.4.7