/home/coin/SVN-release/OS-2.0.0/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   if(numrows == 0) return;
00998   const double * rowLower = getRowLower();
00999   const double * rowUpper = getRowUpper();
01000   rowsense_ = new char [numrows];
01001   rhs_ = new double [numrows];
01002   rowrange_ = new double [numrows];
01003   for(int i = 0 ; i < numrows ; i++) {
01004     rowrange_[i]=0.;
01005     convertBoundToSense(rowLower[i], rowUpper[i], rowsense_[i], rhs_[i], rowrange_[i]);
01006   }
01007 }
01008 
01018 const char *
01019 OsiTMINLPInterface::getRowSense() const
01020 {
01021   if(rowsense_==NULL) {
01022     extractSenseRhsAndRange();
01023   }
01024   return rowsense_;
01025 }
01026 
01035 const double *
01036 OsiTMINLPInterface::getRightHandSide() const
01037 {
01038   if(rhs_==NULL) {
01039     extractSenseRhsAndRange();
01040   }
01041   return rhs_;
01042 }
01043 
01052 const double *
01053 OsiTMINLPInterface::getRowRange() const
01054 {
01055   if(rowrange_==NULL) {
01056     extractSenseRhsAndRange();
01057   }
01058   return rowrange_;
01059 }
01060 
01061 const double *
01062 OsiTMINLPInterface::getRowLower() const
01063 {
01064   return problem_->g_l();
01065 }
01066 
01067 const double *
01068 OsiTMINLPInterface::getRowUpper() const
01069 {
01070   return problem_->g_u();
01071 }
01072 
01074 bool
01075 OsiTMINLPInterface::isContinuous(int colNumber) const
01076 {
01077   return (problem_->var_types()[colNumber]==TMINLP::CONTINUOUS);
01078 }
01079 
01081 bool
01082 OsiTMINLPInterface::isBinary(int colNumber) const
01083 {
01084   return (problem_->var_types()[colNumber]==TMINLP::BINARY);
01085 }
01086 
01091 bool
01092 OsiTMINLPInterface::isInteger(int colNumber) const
01093 {
01094   return ((problem_->var_types()[colNumber]==TMINLP::BINARY)||
01095       (problem_->var_types()[colNumber]==TMINLP::INTEGER));
01096 }
01097 
01099 bool
01100 OsiTMINLPInterface::isIntegerNonBinary(int colNumber) const
01101 {
01102   return (problem_->var_types()[colNumber]==TMINLP::INTEGER);
01103 }
01105 bool
01106 OsiTMINLPInterface::isFreeBinary(int colNumber) const
01107 {
01108   return ((problem_->var_types()[colNumber]==TMINLP::BINARY)
01109       &&((getColUpper()[colNumber]-getColLower()[colNumber]) > 1 - 1e-09));
01110 }
01111 
01113 double
01114 OsiTMINLPInterface::getInfinity() const
01115 {
01116   return COIN_DBL_MAX;
01117 }
01118 
01120 const double *
01121 OsiTMINLPInterface::getColSolution() const
01122 {
01123   if(hasBeenOptimized_)
01124     return problem_->x_sol();
01125   else
01126     return problem_->x_init();
01127 }
01128 
01130 const double *
01131 OsiTMINLPInterface::getRowPrice() const
01132 {
01133   if(hasBeenOptimized_)
01134     return problem_->duals_sol();
01135   else
01136     return problem_->duals_init();
01137 }
01138 
01140 const double *
01141 OsiTMINLPInterface::getReducedCost() const
01142 {
01143   (*handler_)<<"WARNING : trying to access reduced cost in Ipopt always retrun 0"<<CoinMessageEol;
01144   if(reducedCosts_==NULL) {
01145     reducedCosts_ = new double [getNumCols()];
01146     CoinFillN(reducedCosts_,getNumCols(),0.);
01147   }
01148   return reducedCosts_;
01149 }
01150 
01153 const double *
01154 OsiTMINLPInterface::getRowActivity() const
01155 {
01156   return problem_->g_sol();
01157 }
01158 
01162 int
01163 OsiTMINLPInterface::getIterationCount() const
01164 {
01165     return app_->IterationCount();
01166 }
01167 
01168 
01171 void
01172 OsiTMINLPInterface::setColLower( int elementIndex, double elementValue )
01173 {
01174   //  if(fabs(problem_->x_l()[elementIndex]-elementValue)>1e-06)
01175   problem_->SetVariableLowerBound(elementIndex,elementValue);
01176   hasBeenOptimized_ = false;
01177 }
01178 
01181 void
01182 OsiTMINLPInterface::setColUpper( int elementIndex, double elementValue )
01183 {
01184   //  if(fabs(problem_->x_u()[elementIndex]-elementValue)>1e-06)
01185   problem_->SetVariableUpperBound(elementIndex,elementValue);
01186   hasBeenOptimized_ = false;
01187 }
01188 
01191 void
01192 OsiTMINLPInterface::setColLower( const double* array )
01193 {
01194   problem_->SetVariablesLowerBounds(problem_->num_variables(),
01195                                   array);
01196   hasBeenOptimized_ = false;
01197 }
01198 
01201 void
01202 OsiTMINLPInterface::setColUpper( const double* array )
01203 {
01204   problem_->SetVariablesUpperBounds(problem_->num_variables(), 
01205                                   array);
01206   hasBeenOptimized_ = false;
01207 }
01208 
01211 void
01212 OsiTMINLPInterface::setRowLower( int elementIndex, double elementValue )
01213 {
01214   throw SimpleError("Not implemented yet but should be if necessary.",
01215       "setRowLower");
01216   hasBeenOptimized_ = false;
01217 }
01218 
01221 void
01222 OsiTMINLPInterface::setRowUpper( int elementIndex, double elementValue )
01223 {
01224   throw SimpleError("Not implemented yet but should be if necessary.",
01225       "setRowUpper");
01226   hasBeenOptimized_ = false;
01227 }
01228 
01230 void
01231 OsiTMINLPInterface::setRowType(int index, char sense, double rightHandSide,
01232     double range)
01233 {
01234   throw SimpleError("Not implemented yet but should be if necessary.",
01235       "setRowType");
01236   hasBeenOptimized_ = false;
01237 }
01238 
01239 
01242 void
01243 OsiTMINLPInterface::setObjSense(double s)
01244 {
01245   throw SimpleError("Can not change objective sense of an Ipopt problem.",
01246       "setObjSense");
01247   hasBeenOptimized_ = false;
01248 }
01249 
01259 void
01260 OsiTMINLPInterface::setColSolution(const double *colsol)
01261 {
01262   problem_->setxInit(getNumCols(), colsol);
01263   hasBeenOptimized_ = false;
01264 }
01265 
01276 void
01277 OsiTMINLPInterface::setRowPrice(const double * rowprice)
01278 {
01279   problem_->setDualsInit(getNumCols()*2 + getNumRows(), rowprice);
01280   hasBeenOptimized_ = false;
01281 }
01282 
01289 CoinWarmStart *
01290 OsiTMINLPInterface::getEmptyWarmStart () const
01291 {return app_->getEmptyWarmStart();}
01292 
01294 CoinWarmStart* 
01295 OsiTMINLPInterface::getWarmStart() const
01296 {
01297   if (exposeWarmStart_) {
01298     return internal_getWarmStart();;
01299   }
01300   else {
01301     return getEmptyWarmStart();
01302   }
01303 }
01306 bool 
01307 OsiTMINLPInterface::setWarmStart(const CoinWarmStart* ws)
01308 {
01309   if (exposeWarmStart_) {
01310     return internal_setWarmStart(ws);
01311   }
01312   else {
01313     return true;
01314   }
01315 }
01317 CoinWarmStart* 
01318 OsiTMINLPInterface::internal_getWarmStart() const
01319 {
01320   if (exposeWarmStart_ && warmstart_) {
01321     return warmstart_->clone();
01322   }
01323   else {
01324     return getEmptyWarmStart();
01325   }
01326 }
01329 bool 
01330 OsiTMINLPInterface::internal_setWarmStart(const CoinWarmStart* ws)
01331 {
01332   delete warmstart_;
01333   warmstart_ = NULL;
01334   hasBeenOptimized_ = false;
01335   if (exposeWarmStart_) {
01336     if (ws == NULL) {
01337       return true;
01338     }
01339   if(app_->warmStartIsValid(ws)) {
01340     warmstart_ = ws->clone();
01341     return true;
01342   }
01343   // See if it is anything else than the CoinWarmStartBasis that all others
01344   // derive from
01345   const CoinWarmStartPrimalDual* pdws =
01346     dynamic_cast<const CoinWarmStartPrimalDual*>(ws);
01347   if (pdws) {
01348     // Must be an IpoptWarmStart, since the others do not derive from this.
01349     // Accept it.
01350     warmstart_ = new IpoptWarmStart(*pdws);
01351     return true;
01352   }
01353   return false;
01354   }
01355   else {
01356     return true;
01357   }
01358 }
01359 
01361 void
01362 OsiTMINLPInterface::setContinuous(int index)
01363 {
01364   problem_->SetVariableType(index, TMINLP::CONTINUOUS);
01365   hasBeenOptimized_ = false;
01366 }
01368 void
01369 OsiTMINLPInterface::setInteger(int index)
01370 {
01371   problem_->SetVariableType(index, TMINLP::INTEGER);
01372   hasBeenOptimized_ = false;
01373 }
01374 
01376 double
01377 OsiTMINLPInterface::getObjValue() const
01378 {
01379   return problem_->obj_value();
01380 }
01381 
01382 //#############################################################################
01383 // Parameter related methods
01384 //#############################################################################
01385 
01386 bool
01387 OsiTMINLPInterface::setIntParam(OsiIntParam key, int value)
01388 {
01389   //  debugMessage("OsiCpxSolverInterface::setIntParam(%d, %d)\n", key, value);
01390 
01391   bool retval = false;
01392   switch (key) {
01393   case OsiMaxNumIteration:
01394     retval = false;
01395     break;
01396   case OsiMaxNumIterationHotStart:
01397     if( value >= 0 ) {
01398       retval = false;
01399     }
01400     else
01401       retval = false;
01402     break;
01403   case OsiLastIntParam:
01404     retval = false;
01405     break;
01406   default:
01407     retval = false;
01408     (*handler_)<< "Unhandled case in setIntParam\n"<<CoinMessageEol;
01409     break;
01410   }
01411   return retval;
01412 }
01413 
01414 //-----------------------------------------------------------------------------
01415 
01416 bool
01417 OsiTMINLPInterface::setDblParam(OsiDblParam key, double value)
01418 {
01419   //  debugMessage("OsiTMINLPInterface::setDblParam(%d, %g)\n", key, value);
01420 
01421   bool retval = false;
01422   switch (key) {
01423   case OsiDualObjectiveLimit:
01424     OsiDualObjectiveLimit_ = value;
01425     retval = true;
01426     break;
01427   case OsiPrimalObjectiveLimit:
01428     (*handler_)<<"Can not set primal objective limit parameter"<<CoinMessageEol;
01429     retval = false;
01430     break;
01431   case OsiDualTolerance:
01432     (*handler_)<<"Can not set dual tolerance parameter"<<CoinMessageEol;
01433     retval = false;
01434     break;
01435   case OsiPrimalTolerance:
01436     (*handler_)<<"Can not set primal tolerance parameter"<<CoinMessageEol;
01437     retval = false;
01438   case OsiObjOffset:
01439     retval = OsiSolverInterface::setDblParam(key,value);
01440     break;
01441   case OsiLastDblParam:
01442     retval = false;
01443     break;
01444   default:
01445     retval = false;
01446     (*handler_) << "Unhandled case in setDblParam"<<CoinMessageEol;
01447     break;
01448   }
01449   return retval;
01450 }
01451 
01452 
01453 //-----------------------------------------------------------------------------
01454 
01455 bool
01456 OsiTMINLPInterface::setStrParam(OsiStrParam key, const std::string & value)
01457 {
01458   //  debugMessage("OsiTMINLPInterface::setStrParam(%d, %s)\n", key, value.c_str());
01459 
01460   bool retval=false;
01461   switch (key) {
01462   case OsiProbName:
01463     OsiSolverInterface::setStrParam(key,value);
01464     return retval = true;
01465   case OsiSolverName:
01466     return false;
01467   case OsiLastStrParam:
01468     return false;
01469   }
01470   return false;
01471 }
01472 
01473 //-----------------------------------------------------------------------------
01474 
01475 bool
01476 OsiTMINLPInterface::getIntParam(OsiIntParam key, int& value) const
01477 {
01478   //  debugMessage("OsiTMINLPInterface::getIntParam(%d)\n", key);
01479 
01480   value = -COIN_INT_MAX; // Give a dummy value
01481   bool retval = false;
01482   switch (key) {
01483   case OsiMaxNumIteration:
01484     retval = false;
01485     break;
01486   case OsiMaxNumIterationHotStart:
01487     retval = false;
01488     break;
01489   case OsiLastIntParam:
01490     retval = false;
01491     break;
01492   default:
01493     retval = false;
01494     (*handler_) << "Unhandled case in setIntParam"<<CoinMessageEol;
01495   }
01496   return retval;
01497 }
01498 
01499 //-----------------------------------------------------------------------------
01500 
01501 bool
01502 OsiTMINLPInterface::getDblParam(OsiDblParam key, double& value) const
01503 {
01504   //  debugMessage("OsiTMINLPInterface::getDblParam(%d)\n", key);
01505 
01506   bool retval = false;
01507   switch (key) {
01508   case OsiDualObjectiveLimit:
01509     value = OsiDualObjectiveLimit_;
01510     retval = true;
01511     break;
01512   case OsiPrimalObjectiveLimit:
01513     value = getInfinity();
01514     retval = true;
01515     break;
01516   case OsiDualTolerance:
01517     retval = false;
01518     break;
01519   case OsiPrimalTolerance:
01520     options()->GetNumericValue("tol", value,"");
01521     //value = 1e-07;
01522     retval = true;
01523     break;
01524   case OsiObjOffset:
01525     retval = OsiSolverInterface::getDblParam(key, value);
01526     break;
01527   case OsiLastDblParam:
01528     retval = false;
01529     break;
01530   }
01531   return retval;
01532 }
01533 
01534 
01535 //-----------------------------------------------------------------------------
01536 
01537 bool
01538 OsiTMINLPInterface::getStrParam(OsiStrParam key, std::string & value) const
01539 {
01540   //  debugMessage("OsiTMINLPInterface::getStrParam(%d)\n", key);
01541 
01542   switch (key) {
01543   case OsiProbName:
01544     OsiSolverInterface::getStrParam(key, value);
01545     break;
01546   case OsiSolverName:
01547     value = "Ipopt";
01548     break;
01549   case OsiLastStrParam:
01550     return false;
01551   }
01552 
01553   return true;
01554 }
01555 
01556 void
01557 OsiTMINLPInterface::randomStartingPoint()
01558 {
01559   int numcols = getNumCols();
01560   const double * colLower = getColLower();
01561   const double * colUpper = getColUpper();
01562   double * sol = new double[numcols];
01563   const Number * x_init = problem_->x_init_user();
01564   const double* perturb_radius = NULL;
01565   if (randomGenerationType_ == perturb_suffix) {
01566     const TMINLP::PerturbInfo* pertubinfo = tminlp_->perturbInfo();
01567     if (pertubinfo) {
01568       perturb_radius = pertubinfo->GetPerturbationArray();
01569     }
01570     if (!perturb_radius) {
01571       throw SimpleError("Can't use perturb_radius if no radii are given.",
01572                         "randomStartingPoint");
01573     }
01574   }
01575   for(int i = 0 ; i < numcols ; i++) {
01576     int randomGenerationType = randomGenerationType_;
01577     if(x_init[i] < colLower[i] || x_init[i] > colUpper[i])
01578       randomGenerationType = uniform;
01579     if(randomGenerationType_ == uniform){
01580       double lower = std::min(-maxRandomRadius_,colUpper[i] - maxRandomRadius_);
01581       lower = std::max(colLower[i], lower);
01582       double upper = std::max(maxRandomRadius_,colLower[i] + maxRandomRadius_);
01583       upper = std::min(colUpper[i],upper);
01584       lower = std::min(upper,lower);
01585       upper = std::max(upper, lower);
01586       double interval = upper - lower;
01587       sol[i] = CoinDrand48()*(interval) + lower;}
01588     else if (randomGenerationType_ == perturb){
01589       const double lower = std::max(x_init[i] - max_perturbation_, colLower[i]);
01590       const double upper = std::min(x_init[i] + max_perturbation_, colUpper[i]);
01591       const double interval = upper - lower;
01592       sol[i]  = lower + CoinDrand48()*(interval);
01593     }
01594     else if (randomGenerationType_ == perturb_suffix){
01595       const double radius = perturb_radius[i];
01596       const double lower = std::max(x_init[i] - radius*max_perturbation_, colLower[i]);
01597       const double upper = std::min(x_init[i] + radius*max_perturbation_, colUpper[i]);
01598       const double interval = upper - lower;
01599       sol[i]  = lower + CoinDrand48()*(interval);
01600     }
01601   }
01602   app_->disableWarmStart();
01603   setColSolution(sol);
01604   delete [] sol;
01605 }
01606 
01607 
01608 
01610 int OsiTMINLPInterface::initializeJacobianArrays()
01611 {
01612   Index n, m, nnz_h_lag;
01613   TNLP::IndexStyleEnum index_style;
01614   tminlp_->get_nlp_info( n, m, nnz_jac, nnz_h_lag, index_style);
01615 
01616   if(jRow_ != NULL) delete jRow_;
01617   if(jCol_ != NULL) delete jCol_;
01618   if(jValues_ != NULL) delete jValues_;
01619 
01620   jRow_ = new Index[nnz_jac];
01621   jCol_ = new Index[nnz_jac];
01622   jValues_ = new Number[nnz_jac];
01623   tminlp_->eval_jac_g(n, NULL, 0, m, nnz_jac, jRow_, jCol_, NULL);
01624   if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
01625   {
01626     for(int i = 0 ; i < nnz_jac ; i++){
01627       jRow_[i]--;
01628       jCol_[i]--;
01629     }
01630   }
01631 
01632   if(constTypes_ != NULL) delete [] constTypes_;
01633 //  if(constTypesNum_ != NULL) delete [] constTypesNum_;
01634 
01635   constTypes_ = new TNLP::LinearityType[getNumRows()];
01636   tminlp_->get_constraints_linearity(getNumRows(), constTypes_);
01637 //  constTypesNum_ = new int[getNumRows()];
01638   for(int i = 0; i < getNumRows() ; i++) {
01639     if(constTypes_[i]==TNLP::NON_LINEAR) {
01640       //constTypesNum_[i] = 
01641       nNonLinear_++;
01642     }
01643   }
01644   return nnz_jac;
01645 }
01646 
01647 
01648 double 
01649 OsiTMINLPInterface::getConstraintsViolation(const double *x, double &obj)
01650 {
01651   int numcols = getNumCols();
01652   int numrows = getNumRows();
01653   double * g = new double[numrows];
01654   tminlp_->eval_g(numcols, x, 1, numrows, g);
01655   const double * rowLower = getRowLower();
01656   const double * rowUpper = getRowUpper();
01657 
01658 
01659   double norm = 0;
01660   for(int i = 0; i< numrows ; i++) {
01661     if(!constTypes_ || constTypes_[i] == TNLP::NON_LINEAR) {
01662       double rowViolation = 0;
01663       if(rowLower[i] > -1e10)
01664          rowViolation = std::max(0.,rowLower[i] - g[i]);
01665 
01666       if(rowUpper[i] < 1e10);
01667         rowViolation = std::max(rowViolation, g[i] - rowUpper[i]);
01668 
01669       norm = rowViolation > norm ? rowViolation : norm;
01670     }
01671   }
01672   tminlp_->eval_f(numcols, x, 1, obj);
01673   delete [] g;
01674   return norm;
01675 }
01676 
01678 double
01679 OsiTMINLPInterface::getNonLinearitiesViolation(const double *x, const double obj)
01680 {
01681   double f;
01682   double norm = getConstraintsViolation(x, f);
01683   assert((f - obj) > -1e-08);
01684   norm =  (f - obj) > norm ? f - obj : norm;
01685   return norm;
01686 }
01687 
01688 
01689 
01690 //A procedure to try to remove small coefficients in OA cuts (or make it non small
01691 static inline
01692 bool cleanNnz(double &value, double colLower, double colUpper,
01693     double rowLower, double rowUpper, double colsol,
01694     double & lb, double &ub, double tiny, double veryTiny)
01695 {
01696   if(fabs(value)>= tiny) return 1;
01697 
01698   if(fabs(value)<veryTiny) return 0;//Take the risk?
01699 
01700   //try and remove
01701   double infty = 1e20;
01702   bool colUpBounded = colUpper < 10000;
01703   bool colLoBounded = colLower > -10000;
01704   bool rowNotLoBounded =  rowLower <= - infty;
01705   bool rowNotUpBounded = rowUpper >= infty;
01706   bool pos =  value > 0;
01707 
01708   if(colLoBounded && pos && rowNotUpBounded) {
01709     lb += value * (colsol - colLower);
01710     return 0;
01711   }
01712   else
01713     if(colLoBounded && !pos && rowNotLoBounded) {
01714       ub += value * (colsol - colLower);
01715       return 0;
01716     }
01717     else
01718       if(colUpBounded && !pos && rowNotUpBounded) {
01719         lb += value * (colsol - colUpper);
01720         return 0;
01721       }
01722       else
01723         if(colUpBounded && pos && rowNotLoBounded) {
01724           ub += value * (colsol - colUpper);
01725           return 0;
01726         }
01727   //can not remove coefficient increase it to smallest non zero
01728   if(pos) value = tiny;
01729   else
01730     value = - tiny;
01731   return 1;
01732 }
01733 
01736 void
01737 OsiTMINLPInterface::getOuterApproximation(OsiCuts &cs, const double * x, bool getObj, const double * x2, double theta, bool global)
01738 {
01739   int n,m, nnz_jac_g, nnz_h_lag;
01740   TNLP::IndexStyleEnum index_style;
01741   problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
01742   if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
01743     initializeJacobianArrays();
01744   assert(jRow_ != NULL);
01745   assert(jCol_ != NULL);
01746   double * g = new double[m];
01747   problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
01748   problem_to_optimize_->eval_g(n,x,1,m,g);
01749   //As jacobian is stored by cols fill OsiCuts with cuts
01750   CoinPackedVector * cuts = new CoinPackedVector[nNonLinear_ + 1];
01751   double * lb = new double[nNonLinear_ + 1];
01752   double * ub = new double[nNonLinear_ + 1];
01753 
01754   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.
01755   int numCuts = 0;
01756 
01757   const double * rowLower = getRowLower();
01758   const double * rowUpper = getRowUpper();
01759   const double * colLower = getColLower();
01760   const double * colUpper = getColUpper();
01761   const double * duals = getRowPrice() + 2 * n;
01762   double infty = getInfinity();
01763   double nlp_infty = infty_;
01764   
01765   for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
01766     if(constTypes_[rowIdx] == TNLP::NON_LINEAR) {
01767 #if 0
01768       if(fabs(duals[rowIdx]) == 0.)
01769       {
01770         row2cutIdx[rowIdx] = -1;
01771 #ifdef NDEBUG
01772         (*handler_)<<"non binding constraint"<<CoinMessageEol;
01773 #endif
01774         continue;
01775       }
01776 #endif
01777       row2cutIdx[rowIdx] = numCuts;
01778       if(rowLower[rowIdx] > - nlp_infty)
01779         lb[numCuts] = rowLower[rowIdx] - g[rowIdx];
01780       else
01781         lb[numCuts] = - infty;
01782       if(rowUpper[rowIdx] < nlp_infty)
01783         ub[numCuts] = rowUpper[rowIdx] - g[rowIdx];
01784       else
01785         ub[numCuts] = infty;
01786       if(rowLower[rowIdx] > -infty && rowUpper[rowIdx] < infty)
01787       {
01788         if(duals[rowIdx] >= 0)// <= inequality
01789           lb[numCuts] = - infty;
01790         if(duals[rowIdx] <= 0)// >= inequality
01791           ub[numCuts] = infty;
01792       }
01793       
01794       numCuts++;
01795     }
01796     else
01797       row2cutIdx[rowIdx] = -1;
01798   }
01799 
01800 
01801   for(int i = 0 ; i < nnz_jac_g ; i++) {
01802     const int &rowIdx = jRow_[i];
01803     const int & cutIdx = row2cutIdx[ rowIdx ];
01804     if(cutIdx != -1) {
01805       const int &colIdx = jCol_[i];
01806       //"clean" coefficient
01807       if(cleanNnz(jValues_[i],colLower[colIdx], colUpper[colIdx],
01808                   rowLower[rowIdx], rowUpper[rowIdx],
01809                   x[colIdx],
01810                   lb[cutIdx],
01811                   ub[cutIdx], tiny_, veryTiny_)) {
01812         cuts[cutIdx].insert(colIdx,jValues_[i]);
01813         if(lb[cutIdx] > - infty)
01814           lb[cutIdx] += jValues_[i] * x[colIdx];
01815         if(ub[cutIdx] < infty)
01816           ub[cutIdx] += jValues_[i] * x[colIdx];
01817       }
01818     }
01819   }
01820 
01821   int * cut2rowIdx = NULL;
01822   if (IsValid(cutStrengthener_) || oaHandler_->logLevel() > 0) {
01823     cut2rowIdx = new int [numCuts];// Store correspondance between indices of cut and indices of rows. For each cut
01824     for(int rowIdx = 0 ; rowIdx < m ; rowIdx++){
01825        if(row2cutIdx[rowIdx] >= 0){
01826           cut2rowIdx[row2cutIdx[rowIdx]] = rowIdx;
01827        }
01828     }
01829   }
01830 
01831   for(int cutIdx = 0; cutIdx < numCuts; cutIdx++) {
01832     //Compute cut violation
01833     if(x2 != NULL) {
01834       double rhs = cuts[cutIdx].dotProduct(x2);
01835       double violation = 0.;
01836       violation = std::max(violation, rhs - ub[cutIdx]);
01837       violation = std::max(violation, lb[cutIdx] - rhs);
01838       if(violation < theta) {
01839         if(oaHandler_->logLevel() > 0)
01840           oaHandler_->message(CUT_NOT_VIOLATED_ENOUGH, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
01841         continue;}
01842         if(oaHandler_->logLevel() > 0)
01843           oaHandler_->message(VIOLATED_OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
01844     }
01845     else if (oaHandler_->logLevel() > 0)
01846       oaHandler_->message(OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<CoinMessageEol;
01847   OsiRowCut newCut;
01848     //    if(lb[i]>-1e20) assert (ub[i]>1e20);
01849 
01850     if (IsValid(cutStrengthener_)) {
01851       const int& rowIdx = cut2rowIdx[cutIdx];
01852       bool retval =
01853         cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
01854                                        GetRawPtr(problem_), rowIdx,
01855                                        cuts[cutIdx], lb[cutIdx], ub[cutIdx], g[rowIdx],
01856                                        rowLower[rowIdx], rowUpper[rowIdx],
01857                                        n, x, infty);
01858       if (!retval) {
01859         (*messageHandler()) << "error in cutStrengthener_->ComputeCuts\n";
01860         //exit(-2);
01861       }
01862     }
01863     if(global) {
01864       newCut.setGloballyValidAsInteger(1);
01865     }
01866     newCut.setEffectiveness(99.99e99);
01867     newCut.setLb(lb[cutIdx]);
01868     newCut.setUb(ub[cutIdx]);
01869     newCut.setRow(cuts[cutIdx]);
01870     //    CftValidator validator;
01871     //    validator(newCut);
01872     if(oaHandler_->logLevel()>2){
01873       oaHandler_->print(newCut);}
01874     cs.insert(newCut);
01875   }
01876 
01877   delete[] g;
01878   delete [] cuts;
01879   delete [] row2cutIdx;
01880   delete [] cut2rowIdx;
01881 
01882   if(getObj && ! problem_->hasLinearObjective()) { // Get the objective cuts
01883     double * obj = new double [n];
01884     problem_to_optimize_->eval_grad_f(n, x, 1,obj);
01885     double f;
01886     problem_to_optimize_->eval_f(n, x, 1, f);
01887 
01888     CoinPackedVector v;
01889     v.reserve(n);
01890     lb[nNonLinear_] = -f;
01891     ub[nNonLinear_] = -f;
01892     //double minCoeff = 1e50;
01893     for(int i = 0; i<n ; i++) {
01894       if(cleanNnz(obj[i],colLower[i], colUpper[i],
01895           -getInfinity(), 0,
01896           x[i],
01897           lb[nNonLinear_],
01898           ub[nNonLinear_],tiny_, 1e-15)) {
01899         //            minCoeff = std::min(fabs(obj[i]), minCoeff);
01900         v.insert(i,obj[i]);
01901         lb[nNonLinear_] += obj[i] * x[i];
01902         ub[nNonLinear_] += obj[i] * x[i];
01903       }
01904     }
01905     v.insert(n,-1);
01906     //Compute cut violation
01907     bool genCut = true;
01908     if(x2 != NULL) {
01909       double rhs = v.dotProduct(x2);
01910       double violation = std::max(0., rhs - ub[nNonLinear_]);
01911       if(violation < theta) genCut = false;
01912     }
01913     if(genCut) {
01914       if (IsValid(cutStrengthener_)) {
01915         lb[nNonLinear_] = -infty;
01916         bool retval =
01917           cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
01918                                          GetRawPtr(problem_), -1,
01919                                          v, lb[nNonLinear_], ub[nNonLinear_],
01920                                          ub[nNonLinear_], -infty, 0.,
01921                                          n, x, infty);
01922         if (!retval) {
01923     (*handler_)<< "error in cutStrengthener_->ComputeCuts"<<CoinMessageEol;
01924           //exit(-2);
01925         }
01926       }
01927       OsiRowCut newCut;
01928       if(global)
01929         newCut.setGloballyValidAsInteger(1);
01930       newCut.setEffectiveness(99.99e99);
01931       newCut.setRow(v);
01932       newCut.setLb(-COIN_DBL_MAX/*Infinity*/);
01933       newCut.setUb(ub[nNonLinear_]);
01934       //     CftValidator validator;
01935       //     validator(newCut);
01936       cs.insert(newCut);
01937     }
01938     delete [] obj;
01939     }
01940 
01941   delete []lb;
01942   delete[]ub;
01943 }
01944 
01945 
01946 
01949 void
01950 OsiTMINLPInterface::getConstraintOuterApproximation(OsiCuts &cs, int rowIdx, 
01951                                                     const double * x, 
01952                                                     const double * x2, bool global)
01953 {
01954   double g;
01955   int * indices = new int[getNumCols()];
01956   double * values = new double[getNumCols()];
01957   int nnz;
01958   problem_->eval_grad_gi(getNumCols(), x, 1, rowIdx, nnz, indices, values);
01959   problem_->eval_gi(getNumCols(),x,1, rowIdx, g);
01960 
01961   CoinPackedVector cut;
01962   double lb;
01963   double ub;
01964 
01965 
01966   const double rowLower = getRowLower()[rowIdx];
01967   const double rowUpper = getRowUpper()[rowIdx];
01968   const double * colLower = getColLower();
01969   const double * colUpper = getColUpper();
01970   const double dual = (getRowPrice() + 2 * getNumCols())[rowIdx];
01971   double infty = getInfinity();
01972   double nlp_infty = infty_;
01973   
01974   if(rowLower > - nlp_infty)
01975     lb = rowLower - g;
01976   else
01977     lb = - infty;
01978   if(rowUpper < nlp_infty)
01979     ub = rowUpper - g;
01980   else
01981     ub = infty;
01982   if(rowLower > -infty && rowUpper < infty)
01983   {
01984     if(dual >= 0)// <= inequality
01985       lb = - infty;
01986     if(dual <= 0)// >= inequality
01987       ub = infty;
01988   }
01989 
01990   for(int i = 0 ; i < nnz; i++) {
01991      const int &colIdx = indices[i];
01992       //"clean" coefficient
01993       if(cleanNnz(values[i],colLower[colIdx], colUpper[colIdx],
01994                   rowLower, rowUpper,
01995                   x[colIdx],
01996                   lb,
01997                   ub, tiny_, veryTiny_)) {
01998         cut.insert(colIdx,values[i]);
01999         if(lb > - infty)
02000           lb += values[i] * x[colIdx];
02001         if(ub < infty)
02002           ub += values[i] * x[colIdx];
02003     }
02004   }
02005 
02006   OsiRowCut newCut;
02007 
02008   if(global) {
02009     newCut.setGloballyValidAsInteger(1);
02010   }
02011   newCut.setEffectiveness(99.99e99);
02012   newCut.setLb(lb);
02013   newCut.setUb(ub);
02014   newCut.setRow(cut);
02015   cs.insert(newCut);
02016 
02017   delete [] indices;
02018   delete [] values;
02019 }
02020 
02021 void
02022 OsiTMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
02023                                             double a, double s, int L){
02024   if(! IsValid(feasibilityProblem_)) {
02025     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02026   }
02027   feasibilityProblem_->set_use_feasibility_pump_objective(true);
02028   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02029   feasibilityProblem_->setLambda(a);
02030   feasibilityProblem_->setSigma(s);
02031   feasibilityProblem_->setNorm(L);
02032   feasibilityProblem_->set_use_cutoff_constraint(false);
02033   feasibilityProblem_->set_use_local_branching_constraint(false);  
02034   problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
02035   feasibility_mode_ = true;
02036 }
02037 
02038 void
02039 OsiTMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
02040                                                double rhs_local_branching_constraint){
02041   if(! IsValid(feasibilityProblem_)) {
02042     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02043   }
02044   feasibilityProblem_->set_use_feasibility_pump_objective(false);
02045   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02046   feasibilityProblem_->set_use_cutoff_constraint(false);
02047   feasibilityProblem_->set_use_local_branching_constraint(true);  
02048   feasibilityProblem_->set_rhs_local_branching_constraint(rhs_local_branching_constraint);  
02049   problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
02050   feasibility_mode_ = true;
02051 }
02052 
02053 void
02054 OsiTMINLPInterface::switchToOriginalProblem(){
02055   problem_to_optimize_ = GetRawPtr(problem_);
02056   feasibility_mode_ = false;
02057 }
02058 
02059 double
02060 OsiTMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
02061                                             double a, double s, int L)
02062 {
02063   if(! IsValid(feasibilityProblem_)) {
02064     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02065   }
02066   feasibilityProblem_->set_use_feasibility_pump_objective(true);
02067   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02068   feasibilityProblem_->setLambda(a);
02069   feasibilityProblem_->setSigma(s);
02070   feasibilityProblem_->setNorm(L);
02071   feasibilityProblem_->set_use_cutoff_constraint(false);
02072   feasibilityProblem_->set_use_local_branching_constraint(false);  
02073   nCallOptimizeTNLP_++;
02074   totalNlpSolveTime_-=CoinCpuTime();
02075   SmartPtr<TNLPSolver> app2 = app_->clone();
02076   app2->options()->SetIntegerValue("print_level", (Index) 0);
02077   optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
02078   totalNlpSolveTime_+=CoinCpuTime();
02079   hasBeenOptimized_=true;
02080   return getObjValue();
02081 }
02082 
02083 double
02084 OsiTMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
02085                                             int L, double cutoff)
02086 {
02087   if(! IsValid(feasibilityProblem_)) {
02088     throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02089   }
02090   feasibilityProblem_->set_use_feasibility_pump_objective(true);
02091   feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
02092   feasibilityProblem_->setLambda(1.0);
02093   feasibilityProblem_->setSigma(0.0);
02094   feasibilityProblem_->setNorm(L);
02095   feasibilityProblem_->set_use_cutoff_constraint(true);
02096   feasibilityProblem_->set_cutoff(cutoff);
02097   feasibilityProblem_->set_use_local_branching_constraint(false);  
02098   nCallOptimizeTNLP_++;
02099   totalNlpSolveTime_-=CoinCpuTime();
02100   SmartPtr<TNLPSolver> app2 = app_->clone();
02101   app2->options()->SetIntegerValue("print_level", (Index) 0);
02102   optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
02103   totalNlpSolveTime_+=CoinCpuTime();
02104   hasBeenOptimized_=true;
02105   return getObjValue();
02106 }
02107 
02108 double
02109 OsiTMINLPInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs, bool addOnlyViolated, bool global)
02110 {
02111   double ret_val = solveFeasibilityProblem(n, x_bar, inds, 1, 0, 2);
02112   getOuterApproximation(cs, getColSolution(), 0, (addOnlyViolated? x_bar:NULL)
02113                         , global);
02114   return ret_val;
02115 }
02116 
02117 
02118 static bool WarnedForNonConvexOa=false;
02119 
02120 void
02121 OsiTMINLPInterface::extractLinearRelaxation(OsiSolverInterface &si, 
02122                                             const double * x, bool getObj)
02123 {
02124   double * rowLow = NULL;
02125   double * rowUp = NULL;
02126 
02127   int n;
02128   int m;
02129   int nnz_jac_g;
02130   int nnz_h_lag;
02131   TNLP::IndexStyleEnum index_style;
02132   //Get problem information
02133   problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
02134 
02135   //if not allocated allocate spaced for stroring jacobian
02136   if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
02137     initializeJacobianArrays();
02138 
02139   //get Jacobian
02140   problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
02141 
02142 
02143   double *g = new double[m];
02144   problem_to_optimize_->eval_g(n, x, 1, m, g);
02145 
02146   rowLow = new double[m];
02147   rowUp = new double[m];
02148   int * nonBindings = new int[m];//store non binding constraints (which are to be removed from OA)
02149   int numNonBindings = 0;
02150   const double * rowLower = getRowLower();
02151   const double * rowUpper = getRowUpper();
02152   const double * colLower = getColLower();
02153   const double * colUpper = getColUpper();
02154   const double * duals = getRowPrice() + 2*n;
02155   assert(m==getNumRows());
02156   double infty = si.getInfinity();
02157   double nlp_infty = infty_;
02158   
02159   for(int i = 0 ; i < m ; i++) {
02160     if(constTypes_[i] == TNLP::NON_LINEAR) {
02161       //If constraint is range not binding prepare to remove it
02162       if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty && fabs(duals[i]) == 0.)
02163       {
02164         nonBindings[numNonBindings++] = i;
02165         continue;
02166       }
02167       else
02168         if(rowLower[i] > - nlp_infty){
02169           rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
02170           if(! WarnedForNonConvexOa && rowUpper[i] < nlp_infty){
02171              messageHandler()->message(WARNING_NON_CONVEX_OA, messages_)<<CoinMessageEol;
02172              WarnedForNonConvexOa = true;
02173           }
02174         }
02175       else
02176         rowLow[i] = - infty;
02177       if(rowUpper[i] < nlp_infty)
02178         rowUp[i] =  (rowUpper[i] - g[i]) + 1e-07;
02179       else
02180         rowUp[i] = infty;
02181       
02182       //If equality or ranged constraint only add one side by looking at sign of dual multiplier
02183       if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty)
02184       {
02185         if(duals[i] >= 0.)// <= inequality
02186           rowLow[i] = - infty;
02187         if(duals[i] <= 0.)// >= inequality
02188           rowUp[i] = infty;
02189       }
02190     }
02191     else {
02192       if(rowLower[i] > -nlp_infty){
02193       //   printf("Lower %g ", rowLower[i]);
02194          rowLow[i] = (rowLower[i] - g[i]);
02195       }
02196       else
02197         rowLow[i] = - infty;
02198       if(rowUpper[i] < nlp_infty){
02199       //   printf("Upper %g ", rowUpper[i]);
02200          rowUp[i] =  (rowUpper[i] - g[i]);
02201       }
02202       else
02203         rowUp[i] = infty;
02204     }
02205   }
02206 
02207   
02208   
02209   //Then convert everything to a CoinPackedMatrix
02210   //Go through values, clean coefficients and fix bounds
02211   for(int i = 0 ; i < nnz_jac_g ; i++) {
02212     if(constTypes_[jRow_[i]] != TNLP::LINEAR){//For linear just copy is fine.
02213        if(//For other clean tinys
02214        cleanNnz(jValues_[i],colLower[jCol_[i]], colUpper[jCol_[i]],
02215                 rowLower[jRow_[i]], rowUpper[jRow_[i]],
02216                 x[jCol_[i]],
02217                 rowLow[jRow_[i]],
02218                 rowUp[jRow_[i]], tiny_, veryTiny_)) {      
02219           rowLow[jRow_[i]] += jValues_[i] * x[jCol_ [i]];
02220           rowUp[jRow_[i]] += jValues_[i] *x[jCol_[i]];
02221        }
02222     }
02223     else {
02224       double value = jValues_[i] * getColSolution()[jCol_[i]];
02225       rowLow[jRow_[i]] += value;
02226       rowUp[jRow_[i]] += value;
02227     } 
02228   }
02229   CoinPackedMatrix mat(true, jRow_, jCol_, jValues_, nnz_jac_g);
02230   mat.setDimensions(m,n); // In case matrix was empty, this should be enough
02231   
02232   //remove non-bindings equality constraints
02233   mat.deleteRows(numNonBindings, nonBindings);
02234 
02235   int numcols=getNumCols();
02236   double *obj = new double[numcols];
02237   for(int i = 0 ; i < numcols ; i++)
02238     obj[i] = 0.;
02239   
02240   
02241   si.loadProblem(mat, getColLower(), getColUpper(), obj, rowLow, rowUp);
02242   delete [] rowLow;
02243   delete [] rowUp;
02244   delete [] nonBindings;
02245   delete [] g;
02246   for(int i = 0 ; i < getNumCols() ; i++) {
02247     if(isInteger(i))
02248       si.setInteger(i);
02249   }
02250   if(getObj) {
02251      bool addObjVar = false;
02252      if(problem_->hasLinearObjective()){
02253        //Might be in trouble if objective has a constant part
02254        // for now just check that f(0) = 0. 
02255        // If it is not adding a constant term does not seem supported by Osi
02256        // for now
02257        double zero;
02258        problem_to_optimize_->eval_f(n, obj, 1, zero);
02259        si.setDblParam(OsiObjOffset, -zero);
02260        //if(fabs(zero - 0) > 1e-10)
02261          //addObjVar = true;
02262        //else { 
02263          //Copy the linear objective and don't create a dummy variable.
02264          problem_to_optimize_->eval_grad_f(n, x, 1,obj);
02265          si.setObjective(obj);
02266        //}
02267     }
02268     else {
02269       addObjVar = true;
02270    }
02271 
02272    if(addObjVar){
02273       //add variable alpha
02274       //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
02275       CoinPackedVector a;
02276       si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
02277   
02278       // Now get the objective cuts
02279       // get the gradient, pack it and add the cut
02280       problem_to_optimize_->eval_grad_f(n, x, 1,obj);
02281       double ub;
02282       problem_to_optimize_->eval_f(n, x, 1, ub);
02283       ub*=-1;
02284       double lb = -1e300;
02285       CoinPackedVector objCut;
02286       CoinPackedVector * v = &objCut;
02287       v->reserve(n);
02288       for(int i = 0; i<n ; i++) {
02289        if(nnz_jac_g)
02290        {
02291         if(cleanNnz(obj[i],colLower[i], colUpper[i],
02292             -getInfinity(), 0,
02293             x[i],
02294             lb,
02295             ub, tiny_, veryTiny_)) {
02296           v->insert(i,obj[i]);
02297           lb += obj[i] * x[i];
02298           ub += obj[i] * x[i];
02299         }
02300        }
02301        else //Unconstrained problem can not put clean coefficient
02302        {
02303            if(cleanNnz(obj[i],colLower[i], colUpper[i],
02304             -getInfinity(), 0,
02305             x[i],
02306             lb,
02307             ub, 1e-03, 1e-08)) {
02308           v->insert(i,obj[i]);
02309           lb += obj[i] * x[i];
02310           ub += obj[i] * x[i];
02311            }
02312        }
02313       }
02314     v->insert(n,-1);
02315     si.addRow(objCut, lb, ub);
02316     }
02317   }
02318 //  si.writeMpsNative("OA.mps",NULL, NULL, 1);
02319   delete [] obj;
02320 }
02321 
02323 void 
02324 OsiTMINLPInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
02325 { 
02326   if(numberCuts)
02327     freeCachedRowRim();
02328   const OsiRowCut ** cutsPtrs = new const OsiRowCut*[numberCuts];
02329   for(int i = 0 ; i < numberCuts ; i++)
02330   {
02331     cutsPtrs[i] = &cuts[i];
02332   }
02333   problem_->addCuts(numberCuts, cutsPtrs);
02334   delete [] cutsPtrs;
02335 }
02336 
02337 void
02338 OsiTMINLPInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
02339     const char * whereFrom)
02340 {
02341   totalNlpSolveTime_-=CoinCpuTime();
02342   if(warmStarted)
02343     optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
02344   else
02345     optimizationStatus_ = app_->OptimizeTNLP(GetRawPtr(problem_to_optimize_));
02346   totalNlpSolveTime_+=CoinCpuTime();
02347   nCallOptimizeTNLP_++;
02348   hasBeenOptimized_ = true;
02349  
02350  
02351   //Options should have been printed if not done already turn off Ipopt output
02352   if(!hasPrintedOptions) {
02353     hasPrintedOptions = 1;
02354     //app_->Options()->SetIntegerValue("print_level",0, true, true);
02355     app_->options()->SetStringValue("print_user_options","no", false, true);
02356   }
02357   
02358   bool otherDisagree = false ;
02359 #if 0
02360   if(optimizationStatus_ == TNLPSolver::notEnoughFreedom)//Too few degrees of freedom
02361   {
02362     (*messageHandler())<<"Too few degrees of freedom...."<<CoinMessageEol;
02363     int numrows = getNumRows();
02364     int numcols = getNumCols();
02365     
02366     const double * colLower = getColLower();
02367     const double * colUpper = getColUpper();
02368     
02369     
02370     const double * rowLower = getRowLower();
02371     const double * rowUpper = getRowUpper();
02372     
02373     int numberFixed = 0;
02374     for(int i = 0 ; i < numcols ; i++)
02375     {
02376       if(colUpper[i] - colLower[i] <= INT_BIAS)
02377             {
02378               numberFixed++;
02379             }
02380     }
02381     int numberEqualities = 0;
02382     for(int i = 0 ; i < numrows ; i++)
02383     {
02384       if(rowUpper[i] - rowLower[i] <= INT_BIAS)
02385             {
02386               numberEqualities++;
02387             }     
02388     }
02389     if(numcols - numberFixed > numberEqualities || numcols < numberEqualities)
02390     {
02391       std::string probName;
02392       getStrParam(OsiProbName, probName);
02393       throw newUnsolvedError(app_->errorCode(), problem_, probName);
02394     }
02395     double * saveColLow = CoinCopyOfArray(getColLower(), getNumCols());
02396     double * saveColUp = CoinCopyOfArray(getColUpper(), getNumCols());
02397     
02398     for(int i = 0; i < numcols && numcols - numberFixed <= numberEqualities ; i++)
02399     {
02400       if(colUpper[i] - colLower[i] <= INT_BIAS)
02401             {
02402               setColLower(i, saveColLow[i]-1e-06);
02403               setColUpper(i, saveColLow[i]+1e-06);
02404               numberFixed--;
02405             }
02406     }
02407     solveAndCheckErrors(warmStarted, throwOnFailure, whereFrom);
02408     //restore
02409     for(int i = 0; i < numcols && numcols - numberFixed < numrows ; i++)
02410     {
02411       problem_->SetVariableLowerBound(i,saveColLow[i]);
02412       problem_->SetVariableUpperBound(i,saveColUp[i]);
02413     }
02414     delete [] saveColLow;
02415     delete [] saveColUp;
02416     return;
02417   }
02418   else 
02419 #endif
02420     if(!app_->isRecoverable(optimizationStatus_))//Solver failed and the error can not be recovered, throw it
02421     {
02422       std::string probName;
02423       getStrParam(OsiProbName, probName);
02424       throw newUnsolvedError(app_->errorCode(), problem_, probName);
02425     }
02426     else if(testOthers_ && !app_->isError(optimizationStatus_)){
02427       Ipopt::SmartPtr<TMINLP2TNLP> problem_copy = problem_->clone();
02428       //Try other solvers and see if they agree
02429       int f =1;
02430       for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
02431           i != debug_apps_.end() ; i++){
02432         TNLPSolver::ReturnStatus otherStatus = (*i)->OptimizeTNLP(GetRawPtr(problem_copy));
02433        messageHandler()->message(LOG_LINE, messages_)
02434          <<'d'<<f++<<statusAsString(otherStatus)<<problem_copy->obj_value()
02435          <<(*i)->IterationCount()<<(*i)->CPUTime()<<CoinMessageEol;
02436         if(!(*i)->isError(otherStatus)){
02437            CoinRelFltEq eq(1e-05);
02438            if(otherStatus != optimizationStatus_){
02439              otherDisagree = true;
02440              messageHandler()->message(SOLVER_DISAGREE_STATUS, messages_)
02441              <<app_->solverName()<<statusAsString()
02442              <<(*i)->solverName()<<statusAsString(otherStatus)<<CoinMessageEol; 
02443            }
02444            else if(isProvenOptimal() && !eq(problem_->obj_value(),problem_copy->obj_value()))
02445            {
02446              otherDisagree = true;
02447              messageHandler()->message(SOLVER_DISAGREE_VALUE, messages_)
02448              <<app_->solverName()<<problem_->obj_value()
02449              <<(*i)->solverName()<<problem_copy->obj_value()<<CoinMessageEol; 
02450            }
02451         }
02452      }
02453   }
02454   try{
02455     totalIterations_ += app_->IterationCount();
02456   }
02457   catch(SimpleError &E)
02458   {
02459     if (throwOnFailure)//something failed throw
02460     {
02461       throw SimpleError("No statistics available from Ipopt",whereFrom);
02462     }
02463     else {
02464       return;
02465     }
02466   }
02467   if(problem_->hasUpperBoundingObjective()){//Check if solution is integer and recompute objective value using alternative objective function
02468     const double * sol = getColSolution();
02469     bool integerSol = true;
02470     double intTol = 1e-08;
02471     if(objects()){
02472       int nObjects = numberObjects();
02473       OsiObject ** object = objects();
02474       for(int i = 0 ; i< nObjects ; i++){
02475         int dummy;
02476         if(object[i]->infeasibility(this,dummy) > intTol)
02477         {
02478           integerSol=false;
02479           break;
02480         }
02481       }
02482     }
02483     else{//Only works for integer constraints
02484       int numcols = getNumCols();
02485       for(int i = 0 ; i < numcols ; i++){
02486         if(isInteger(i) || isBinary(i)){
02487           if(fabs(sol[i] - floor(sol[i]+0.5)) > intTol){
02488             integerSol = false;
02489             break;
02490           }
02491         }
02492       }
02493     }
02494     if(integerSol&&isProvenOptimal()){
02495       double help= problem_->evaluateUpperBoundingFunction(sol);
02496      
02497 
02498       OsiAuxInfo * auxInfo = getAuxiliaryInfo();
02499       Bonmin::AuxInfo * bonInfo = dynamic_cast<Bonmin::AuxInfo *>(auxInfo);
02500       if(bonInfo!=0)
02501       {
02502         
02503         if(help<bonInfo->bestObj2())
02504         {
02505           bonInfo->setBestObj2(help);
02506           bonInfo->setBestSolution2(getNumCols(), const_cast<double *>(getColSolution()));
02507 
02508            messageHandler()->message(ALTERNATE_OBJECTIVE, messages_)
02509            <<help<<CoinMessageEol;
02510         }
02511       }
02512       else {
02513         printf("\nWARNING: the algorithm selected does not consider the second objective function\n");
02514       }
02515     }
02516   }
02517   messageHandler()->message(IPOPT_SUMMARY, messages_)
02518     <<whereFrom<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
02519   
02520   if((nCallOptimizeTNLP_ % 20) == 1)
02521     messageHandler()->message(LOG_HEAD, messages_)<<CoinMessageEol;
02522   
02523   
02524   if ( (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) ||
02525        ( otherDisagree )){
02526     messageHandler()->message(SUSPECT_PROBLEM,
02527                               messages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
02528     std::string subProbName;
02529     getStrParam(OsiProbName, subProbName);
02530     std::ostringstream os;
02531     os<<"_"<<nCallOptimizeTNLP_;
02532     subProbName+=os.str();
02533     problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
02534   }
02535   
02536 }
02537 
02539 // Solve Methods                                                  //
02542 void OsiTMINLPInterface::initialSolve()
02543 {
02544   assert(IsValid(app_));
02545   assert(IsValid(problem_));
02546 
02547   // Discard warmstart_ if we had one
02548   delete warmstart_;
02549   warmstart_ = NULL;
02550   
02551   if(!hasPrintedOptions) {
02552     int printOptions;
02553     app_->options()->GetEnumValue("print_user_options",printOptions,"bonmin.");
02554     if(printOptions)
02555       app_->options()->SetStringValue("print_user_options","yes",true,true);
02556   }
02557   if(exposeWarmStart_)
02558     app_->disableWarmStart(); 
02559   solveAndCheckErrors(0,1,"initialSolve");
02560   
02561   //Options should have been printed if not done already turn off Ipopt output
02562   if(!hasPrintedOptions) {
02563     hasPrintedOptions = 1;
02564     app_->options()->SetStringValue("print_user_options","no");
02565     app_->options()->SetIntegerValue("print_level",0);
02566   }
02567   
02568   messageHandler()->message(LOG_FIRST_LINE, messages_)<<nCallOptimizeTNLP_
02569                                                       <<statusAsString()
02570                                                       <<getObjValue()
02571                                                       <<app_->IterationCount()
02572                                                       <<app_->CPUTime()
02573                                                       <<CoinMessageEol;
02574   
02575   int numRetry = firstSolve_ ? numRetryInitial_ : numRetryResolve_;
02576   if(isAbandoned()) {
02577     resolveForRobustness(numRetryUnsolved_);
02578   }
02579   else if(numRetry)
02580     {
02581       resolveForCost(numRetry, numRetryInitial_ > 0);
02583       numRetryInitial_ = 0;
02584     }
02585   firstSolve_ = false;
02586 
02587   // if warmstart_ is not empty then had to use resolveFor... and that created
02588   // the warmstart at the end, and we have nothing to do here. Otherwise...
02589   if (! warmstart_ && ! isAbandoned()) {
02590     if (exposeWarmStart_) {
02591       warmstart_ = app_->getWarmStart(problem_);
02592     }
02593   }
02594 }
02595 
02598 void
02599 OsiTMINLPInterface::resolve()
02600 {
02601   assert(IsValid(app_));
02602   assert(IsValid(problem_));
02603   
02604   int has_warmstart = warmstart_ == NULL ? 0 : 1;
02605   if(warmstart_ == NULL) has_warmstart = 0;
02606   else if(!app_->warmStartIsValid(warmstart_)) has_warmstart = 1;
02607   else has_warmstart = 2;
02608   if (has_warmstart < 2) {
02609     // empty (or unrecognized) warmstart
02610     initialSolve();
02611     return;
02612   }
02613   app_->setWarmStart(warmstart_, problem_);
02614   delete warmstart_;
02615   warmstart_ = NULL;
02616 
02617   if (INT_BIAS > 0.) {
02618     app_->options()->SetStringValue("warm_start_same_structure", "yes");
02619   }
02620   else {
02621     app_->options()->SetStringValue("warm_start_same_structure", "no");
02622   }
02623 
02624   if(problem_->duals_init() != NULL)
02625     app_->enableWarmStart();
02626   else app_->disableWarmStart();
02627   solveAndCheckErrors(1,1,"resolve");
02628   
02629   messageHandler()->message(LOG_FIRST_LINE, messages_)<<nCallOptimizeTNLP_
02630                                                       <<statusAsString()
02631                                                       <<getObjValue()
02632                                                       <<app_->IterationCount()
02633                                                       <<app_->CPUTime()<<CoinMessageEol;
02634   
02635   if(isAbandoned()) {
02636     resolveForRobustness(numRetryUnsolved_);
02637   }
02638   else if(numRetryResolve_ ||
02639           (numRetryInfeasibles_ && isProvenPrimalInfeasible() ))
02640     resolveForCost(std::max(numRetryResolve_, numRetryInfeasibles_), 0);
02641 
02642   // if warmstart_ is not empty then had to use resolveFor... and that created
02643   // the warmstart at the end, and we have nothing to do here. Otherwise...
02644   if (! warmstart_ && ! isAbandoned()) {
02645     if (exposeWarmStart_) {
02646       warmstart_ = app_->getWarmStart(problem_);
02647     }
02648   }
02649 }
02650 
02651 
02653 // Methods returning info on how the solution process terminated  //
02656 bool OsiTMINLPInterface::isAbandoned() const
02657 {
02658   return (
02659         (optimizationStatus_==TNLPSolver::iterationLimit)||
02660         (optimizationStatus_==TNLPSolver::computationError)||
02661         (optimizationStatus_==TNLPSolver::illDefinedProblem)||
02662         (optimizationStatus_==TNLPSolver::illegalOption)||
02663         (optimizationStatus_==TNLPSolver::externalException)|
02664         (optimizationStatus_==TNLPSolver::exception)
02665       );
02666 }
02667 
02669 bool OsiTMINLPInterface::isProvenOptimal() const
02670 {
02671   return (optimizationStatus_==TNLPSolver::solvedOptimal) ||
02672           (optimizationStatus_==TNLPSolver::solvedOptimalTol);
02673 }
02675 bool OsiTMINLPInterface::isProvenPrimalInfeasible() const
02676 {
02677   return (optimizationStatus_ == TNLPSolver::provenInfeasible);
02678 }
02680 bool OsiTMINLPInterface::isProvenDualInfeasible() const
02681 {
02682   return (optimizationStatus_ == TNLPSolver::unbounded);
02683 }
02685 bool OsiTMINLPInterface::isPrimalObjectiveLimitReached() const
02686 {
02687   (*handler_)<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<CoinMessageEol;
02688   return 0;
02689 }
02691 bool OsiTMINLPInterface::isDualObjectiveLimitReached() const
02692 {
02693   //  (*messageHandler_)<<"Warning : isDualObjectiveLimitReached not implemented yet"<<CoinMessageEol;
02694   return (optimizationStatus_==TNLPSolver::unbounded);
02695 
02696 }
02698 bool OsiTMINLPInterface::isIterationLimitReached() const
02699 {
02700   return (optimizationStatus_==TNLPSolver::iterationLimit);
02701 }
02702 
02703 void
02704 OsiTMINLPInterface::extractInterfaceParams()
02705 {
02706   if (IsValid(app_)) {
02707     int logLevel;
02708     app_->options()->GetIntegerValue("nlp_log_level", logLevel,"bonmin.");
02709     messageHandler()->setLogLevel(logLevel);
02710 
02711 #ifdef COIN_HAS_FILTERSQP
02712     FilterSolver * filter = dynamic_cast<FilterSolver *>(GetRawPtr(app_));
02713 
02714     bool is_given =
02715 #endif
02716       app_->options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,"bonmin.");
02717 
02718 #ifdef COIN_HAS_FILTERSQP
02719     if(filter && !is_given){
02720       // overwriting default for filter
02721       maxRandomRadius_ = 10.;
02722     }
02723 #endif
02724    
02725    int oaCgLogLevel = 0;
02726    app_->options()->GetIntegerValue("oa_cuts_log_level", oaCgLogLevel,"bonmin.");
02727    oaHandler_->setLogLevel(oaCgLogLevel); 
02728     
02729     int exposeWs = false;
02730     app_->options()->GetEnumValue("warm_start", exposeWs, "bonmin.");
02731     setExposeWarmStart(exposeWs > 0);
02732     
02733     app_->options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,"bonmin.");
02734     app_->options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,"bonmin.");
02735     app_->options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,"bonmin.");
02736     app_->options()->GetIntegerValue("num_resolve_at_infeasibles", numRetryInfeasibles_,"bonmin.");
02737     app_->options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,"bonmin.");
02738     app_->options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,"bonmin.");
02739     app_->options()->GetNumericValue
02740     ("warm_start_bound_frac" ,pushValue_,"bonmin.");
02741     app_->options()->GetNumericValue("tiny_element",tiny_,"bonmin.");
02742     app_->options()->GetNumericValue("very_tiny_element",veryTiny_,"bonmin.");
02743     app_->options()->GetNumericValue("random_point_perturbation_interval",max_perturbation_,"bonmin.");
02744     app_->options()->GetEnumValue("random_point_type",randomGenerationType_,"bonmin.");
02745     int cut_strengthening_type;
02746     app_->options()->GetEnumValue("cut_strengthening_type", cut_strengthening_type,"bonmin.");
02747 
02748     if (cut_strengthening_type != CS_None) {
02749       // TNLP solver to be used in the cut strengthener
02750       cutStrengthener_ = new CutStrengthener(app_->clone(), app_->options());
02751     }
02752   }
02753 }
02754 
02755 void
02756 OsiTMINLPInterface::SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)
02757 {
02758   strong_branching_solver_ = strong_branching_solver;
02759 }
02760 
02761   //#define STRONG_COMPARE
02762 #ifdef STRONG_COMPARE
02763   static double objorig;
02764 #endif
02765 
02766 void
02767 OsiTMINLPInterface::markHotStart()
02768 {
02769   if (IsValid(strong_branching_solver_)) {
02770 #ifdef STRONG_COMPARE
02771     // AWDEBUG
02772     OsiSolverInterface::markHotStart();
02773     objorig = getObjValue();
02774 #endif
02775     optimizationStatusBeforeHotStart_ = optimizationStatus_;
02776     strong_branching_solver_->markHotStart(this);
02777   }
02778   else {
02779     // Default Implementation
02780     OsiSolverInterface::markHotStart();
02781   }
02782 }
02783 
02784 void
02785 OsiTMINLPInterface::solveFromHotStart()
02786 {
02787   if (IsValid(strong_branching_solver_)) {
02788 #ifdef STRONG_COMPARE
02789     // AWDEBUG
02790     OsiSolverInterface::solveFromHotStart();
02791     double obj_nlp = getObjValue() - objorig;
02792 #endif
02793     optimizationStatus_ = strong_branching_solver_->solveFromHotStart(this);
02794     hasBeenOptimized_ = true;
02795 #ifdef STRONG_COMPARE
02796     double obj_other = getObjValue() - objorig;
02797     printf("AWDEBUG: Strong Branching results: NLP = %15.8e Other = %15.8e\n",
02798            obj_nlp, obj_other);
02799 #endif
02800   }
02801   else {
02802     // Default Implementation
02803     OsiSolverInterface::solveFromHotStart();
02804   }
02805 }
02806 
02807 void
02808 OsiTMINLPInterface::unmarkHotStart()
02809 {
02810   if (IsValid(strong_branching_solver_)) {
02811 #ifdef STRONG_COMPARE
02812     // AWDEBUG
02813     OsiSolverInterface::unmarkHotStart();
02814 #endif
02815     strong_branching_solver_->unmarkHotStart(this);
02816     optimizationStatus_ = optimizationStatusBeforeHotStart_;
02817   }
02818   else {
02819     // Default Implementation
02820     OsiSolverInterface::unmarkHotStart();
02821   }
02822 }
02823 
02824 const double * OsiTMINLPInterface::getObjCoefficients() const
02825 {
02826   const int n = getNumCols();
02827   delete [] obj_;
02828   obj_ = NULL;
02829   obj_ = new double[n];
02830 
02831   bool new_x = true;
02832   const double* x_sol = problem_->x_sol();
02833   bool retval = problem_->eval_grad_f(n, x_sol, new_x, obj_);
02834   
02835   if (!retval) {
02836     // Let's see if that happens - it will cause a crash
02837     fprintf(stderr, "ERROR WHILE EVALUATING GRAD_F in OsiTMINLPInterface::getObjCoefficients()\n");
02838     delete [] obj_;
02839     obj_ = NULL;
02840   }
02841 
02842   return obj_;
02843 }
02844 
02845 
02846 }

Generated on Mon Aug 3 03:02:18 2009 by  doxygen 1.4.7