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

Generated on Thu Sep 22 03:05:55 2011 by  doxygen 1.4.7