source: trunk/Bonmin/src/Interfaces/BonOsiTMINLPInterface.cpp @ 1512

Last change on this file since 1512 was 1512, checked in by pbonami, 11 years ago

For ticket 38

File size: 91.0 KB
Line 
1// (C) Copyright International Business Machines Corporation and
2// Carnegie Mellon University 2004, 2007
3//
4// All Rights Reserved.
5// This code is published under the Common Public License.
6//
7// Authors :
8// Pierre Bonami, Carnegie Mellon University,
9// Carl D. Laird, Carnegie Mellon University,
10// Andreas Waechter, International Business Machines Corporation
11//
12// Date : 12/01/2004
13
14#include "BonminConfig.h"
15
16#include "BonOsiTMINLPInterface.hpp"
17#include "BonTMINLP.hpp"
18#include "BonTMINLP2TNLP.hpp"
19#include "BonTNLP2FPNLP.hpp"
20#include "BonTNLPSolver.hpp"
21#include "CoinTime.hpp"
22#include <climits>
23#include <string>
24#include <sstream>
25#include "BonAuxInfos.hpp"
26
27#include "Ipopt/BonIpoptSolver.hpp"
28#include "Ipopt/BonIpoptWarmStart.hpp"
29#ifdef COIN_HAS_FILTERSQP
30#include "Filter/BonFilterSolver.hpp"
31#include "Filter/BonFilterWarmStart.hpp"
32//#include "Filter/BonBqpdWarmStart.hpp"
33#endif
34
35#include "OsiBranchingObject.hpp"
36#include "BonStrongBranchingSolver.hpp"
37
38//Macros to try and make messages definition less heavy
39#include "BonMsgUtils.hpp"
40
41using namespace Ipopt;
42
43
44namespace Bonmin {
45///Register options
46static void
47register_general_options
48(SmartPtr<RegisteredOptions> roptions)
49{
50  roptions->SetRegisteringCategory("nlp interface option", RegisteredOptions::BonminCategory);
51  roptions->AddStringOption3("nlp_solver",
52                             "Choice of the solver for local optima of continuous nlp's",
53                             "Ipopt",
54                             "Ipopt", "Interior Point OPTimizer (https://projects.coin-or.org/Ipopt)",
55                             "filterSQP", "Sequential quadratic programming trust region "
56                                          "algorithm (http://www-unix.mcs.anl.gov/~leyffer/solvers.html)",
57                             "all", "run all available solvers at each node",
58                             "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.");
59  roptions->setOptionExtraInfo("nlp_solver",15);
60  roptions->AddBoundedIntegerOption("nlp_log_level",
61                                    "specify NLP solver interface log level (independent from ipopt print_level).",
62                                     0,2,1,
63                                    "Set the level of output of the OsiTMINLPInterface : "
64                                    "0 - none, 1 - normal, 2 - verbose"
65                                   );
66  roptions->setOptionExtraInfo("nlp_log_level",15);
67
68  roptions->AddStringOption2("file_solution",
69       "Write a file bonmin.sol with the solution",
70       "no",
71       "yes","",
72       "no","","");
73
74  roptions->AddStringOption3("warm_start",
75      "Select the warm start method",
76      "none",
77      "none","No warm start",
78      "optimum","Warm start with direct parent optimum",
79      "interior_point","Warm start with an interior point of direct parent",
80      "This will affect the function getWarmStart(), and as a consequence the warm starting in the various algorithms.");
81  roptions->setOptionExtraInfo("warm_start",8);
82
83  roptions->SetRegisteringCategory("Nlp solution robustness", RegisteredOptions::BonminCategory);
84
85  roptions->AddLowerBoundedNumberOption("max_random_point_radius",
86      "Set max value r for coordinate of a random point.",
87      0.,1,1e5,
88      "When picking a random point coordinate i will be in the interval [min(max(l,-r),u-r), max(min(u,r),l+r)] "
89      "(where l is the lower bound for the variable and u is its upper bound)");
90  roptions->setOptionExtraInfo("max_random_point_radius",8);
91
92  roptions->AddStringOption3("random_point_type","method to choose a random starting point",
93                             "Jon",
94                             "Jon", "Choose random point uniformly between the bounds",
95                             "Andreas", "perturb the starting point of the problem within a prescribed interval",
96                             "Claudia", "perturb the starting point using the perturbation radius suffix information",
97                             "");
98  roptions->setOptionExtraInfo("random_point_type",8);
99
100    roptions->AddLowerBoundedNumberOption("random_point_perturbation_interval",
101                                           "Amount by which starting point is perturbed when choosing to pick random point by perturbating starting point",
102                                           0.,true, 1.,
103                                           "");
104  roptions->setOptionExtraInfo("random_point_perturbation_interval",8);
105                                           
106
107  roptions->AddLowerBoundedIntegerOption
108  ("num_iterations_suspect",
109   "Number of iterations over which a node is considered \"suspect\" (for debugging purposes only, see detailed documentation).",
110   -1,-1,
111   "When the number of iterations to solve a node is above this number, the subproblem at this"
112   " node is considered to be suspect and it will be outputed in a file (set to -1 to deactivate this).");
113  roptions->setOptionExtraInfo("num_iterations_suspect",15);
114
115 
116
117  roptions->AddLowerBoundedIntegerOption("num_retry_unsolved_random_point",
118      "Number $k$ of times that the algorithm will try to resolve an unsolved NLP with a random starting point "
119      "(we call unsolved an NLP for which Ipopt is not "
120      "able to guarantee optimality within the specified tolerances).",
121      0,0,
122      "When Ipopt fails to solve a continuous NLP sub-problem, if $k > 0$, the algorithm will "
123      "try again to solve the failed NLP with $k$ new randomly chosen starting points "
124      " or until the problem is solved with success.");
125  roptions->setOptionExtraInfo("num_retry_unsolved_random_point",15);
126
127
128  roptions->SetRegisteringCategory("Options for non-convex problems", RegisteredOptions::BonminCategory);
129
130
131  roptions->AddLowerBoundedIntegerOption("num_resolve_at_root",
132      "Number $k$ of tries to resolve the root node with different starting points.",
133      0,0,
134      "The algorithm will solve the root node with $k$ random starting points"
135      " and will keep the best local optimum found.");
136  roptions->setOptionExtraInfo("num_resolve_at_root",8);
137
138  roptions->AddLowerBoundedIntegerOption("num_resolve_at_node",
139      "Number $k$ of tries to resolve a node (other than the root) of the tree with different starting point.",
140      0,0,
141      "The algorithm will solve all the nodes with $k$ different random starting points "
142      "and will keep the best local optimum found.");
143  roptions->setOptionExtraInfo("num_resolve_at_node",8);
144
145  roptions->AddLowerBoundedIntegerOption("num_resolve_at_infeasibles",
146      "Number $k$ of tries to resolve an infeasible node (other than the root) of the tree with different starting point.",
147      0,0,
148      "The algorithm will solve all the infeasible nodes with $k$ different random starting points "
149      "and will keep the best local optimum found.");
150  roptions->setOptionExtraInfo("num_resolve_at_infeasibles",8);
151
152
153
154  }
155
156static void register_OA_options
157(SmartPtr<RegisteredOptions> roptions)
158{
159  roptions->SetRegisteringCategory("Outer Approximation cuts generation", RegisteredOptions::BonminCategory);
160 
161  roptions->AddStringOption2("disjunctive_cut_type",
162      "Determine if and what kind of disjunctive cuts should be computed.",
163      "none",
164      "none", "No disjunctive cuts.",
165      "most-fractional", "If discrete variables present, compute disjunction for most-fractional variable");
166  roptions->setOptionExtraInfo("disjunctive_cut_type",7);
167
168  roptions->AddStringOption2("oa_cuts_scope","Specify if OA cuts added are to be set globally or locally valid",
169                             "global",
170                             "local","Cuts are treated as globally valid",
171                             "global", "Cuts are treated as locally valid",
172                             "");
173  roptions->setOptionExtraInfo("oa_cuts_scope",7);
174
175  roptions->AddStringOption2("add_only_violated_oa","Do we add all OA cuts or only the ones violated by current point?",
176                             "no",
177                             "no","Add all cuts",
178                             "yes","Add only violated Cuts","");
179  roptions->setOptionExtraInfo("add_only_violated_oa",7);
180
181  roptions->AddStringOption4("cut_strengthening_type",
182                             "Determines if and what kind of cut strengthening should be performed.",
183                             "none",
184                             "none", "No strengthening of cuts.",
185                             "sglobal", "Strengthen global cuts.",
186                             "uglobal-slocal", "Unstrengthened global and strengthened local cuts",
187                             "sglobal-slocal", "Strengthened global and strengthened local cuts",
188                             "");
189  roptions->setOptionExtraInfo("cut_strengthening_type",7);
190 
191  roptions->AddLowerBoundedNumberOption("tiny_element","Value for tiny element in OA cut",
192      -0.,0,1e-08,
193      "We will remove \"cleanly\" (by relaxing cut) an element lower"
194      " than this.");
195  roptions->setOptionExtraInfo("tiny_element",7);
196
197  roptions->AddLowerBoundedNumberOption("very_tiny_element","Value for very tiny element in OA cut",
198      -0.,0,1e-17,
199      "Algorithm will take the risk of neglecting an element lower"
200      " than this.");
201  roptions->setOptionExtraInfo("very_tiny_element",7);
202
203  roptions->AddLowerBoundedIntegerOption("oa_cuts_log_level",
204                                         "level of log when generating OA cuts.",
205                                         0, 0,
206                                         "0: outputs nothing,\n"
207                                         "1: when a cut is generated, its violation and index of row from which it originates,\n"
208                                         "2: always output violation of the cut.\n"
209                                         "3: output generated cuts incidence vectors.");
210  roptions->setOptionExtraInfo("oa_cuts_log_level",7);
211
212}
213
214
215///Register options
216void
217OsiTMINLPInterface::registerOptions
218(SmartPtr<RegisteredOptions> roptions)
219{
220  // We try to register the options - if those have been registered
221  // already, we catch the exception and don't need to do it again
222  try {
223    register_general_options(roptions);
224    register_OA_options(roptions);
225#ifdef COIN_HAS_FILTERSQP
226    FilterSolver::RegisterOptions(roptions);
227#endif
228#ifdef COIN_HAS_IPOPT
229    IpoptSolver::RegisterOptions(roptions);
230#endif
231  }   
232  catch(RegisteredOptions::OPTION_ALREADY_REGISTERED) {
233    // skipping
234  }
235
236
237}
238
239/** To set some application specific defaults. */
240void 
241OsiTMINLPInterface::setAppDefaultOptions(SmartPtr<OptionsList> Options)
242{}
243
244OsiTMINLPInterface::Messages::Messages
245():CoinMessages((int)OSITMINLPINTERFACE_DUMMY_END)
246{
247  strcpy(source_ ,"NLP");
248  ADD_MSG(SOLUTION_FOUND, std_m, 2,
249          "After %d tries found a solution of %g (previous best %g).");
250
251  ADD_MSG(INFEASIBLE_SOLUTION_FOUND, std_m, 2,
252          "After %d tries found an solution of %g infeasible problem.");
253
254  ADD_MSG(UNSOLVED_PROBLEM_FOUND, std_m, 2,
255          "After %d tries found an solution of %g unsolved problem.");
256  ADD_MSG(WARN_SUCCESS_WS, warn_m, 2,
257       "Problem not solved with warm start but solved without");
258
259  ADD_MSG(WARNING_RESOLVING, warn_m,2,
260       "Trying to resolve NLP with different starting point (%d attempts).");
261  ADD_MSG(WARN_SUCCESS_RANDOM, warn_m, 1,
262       "Problem initially not solved but solved with a random starting point (success on %d attempt)");
263  ADD_MSG(WARN_CONTINUING_ON_FAILURE, warn_m, 1,
264       "Warning : continuing branching, while there are unrecovered failures at nodes");
265
266  ADD_MSG(SUSPECT_PROBLEM, std_m, 2, "NLP number %d is suspect (see bounds and start file)");
267  ADD_MSG(IPOPT_SUMMARY, std_m, 2, "Ipopt return (for %s): status %2d, iter count %4d, time %g");
268  ADD_MSG(BETTER_SOL, std_m, 2, "Solution of value %g found on %d'th attempt");
269
270  ADD_MSG(LOG_HEAD, std_m, 1,
271          "\n          "
272          "    Num      Status      Obj             It       time                 Location");
273  ADD_MSG(LOG_LINE, std_m, 1, 
274          "%c    %-8d %-11s %-15.5g %-8d %-8g %20s");
275
276  ADD_MSG(ALTERNATE_OBJECTIVE, std_m, 1, "Objective value recomputed with alternate objective: %g.");
277 
278  ADD_MSG(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, warn_m, 1, 
279         "resolve called before any call to initialSol  can not use warm starts.");
280  ADD_MSG(ERROR_NO_TNLPSOLVER, warn_m, 1,"Can not parse options when no IpApplication has been created");
281  ADD_MSG(WARNING_NON_CONVEX_OA, warn_m, 1,
282          "OA on non-convex constraint is very experimental.");                         
283  ADD_MSG(SOLVER_DISAGREE_STATUS, warn_m, 1, "%s says problem %s, %s says %s.");
284  ADD_MSG(SOLVER_DISAGREE_VALUE, warn_m, 1, "%s gives objective %.16g, %s gives %.16g.");
285
286}
287
288
289void 
290OsiTMINLPInterface::OaMessageHandler::print(OsiRowCut &row){
291  FILE * fp = filePointer();
292  const int &n = row.row().getNumElements();
293  fprintf(fp,"Row cut has %d elements. Lower bound: %g, upper bound %g.\n", n, row.lb(), row.ub());
294  const int * idx = row.row().getIndices();
295  const double * val = row.row().getElements();
296  for(int i = 0 ; i < n ; i++){
297    fprintf(fp,"%g, x%d",val[i], idx[i]);
298    if(i && i % 7 == 0)
299      fprintf(fp,"\n");
300  } 
301}
302
303OsiTMINLPInterface::OaMessages::OaMessages(): CoinMessages((int) OA_MESSAGES_DUMMY_END){
304   strcpy(source_,"OaCg");
305   ADD_MSG(VIOLATED_OA_CUT_GENERATED, std_m, 1,"Row %d, cut violation is %g: Outer approximation cut generated.");
306
307   ADD_MSG(CUT_NOT_VIOLATED_ENOUGH, std_m, 2,"Row %d, cut violation is %g: Outer approximation cut not generated.");
308
309   ADD_MSG(OA_CUT_GENERATED, std_m, 1,"Row %d: Outer approximation cut not generated.");
310}
311bool OsiTMINLPInterface::hasPrintedOptions=0;
312
313////////////////////////////////////////////////////////////////////
314// Constructors and desctructors                                  //
315////////////////////////////////////////////////////////////////////
316/// Default Constructor
317OsiTMINLPInterface::OsiTMINLPInterface():
318    OsiSolverInterface(),
319    tminlp_(NULL),
320    problem_(NULL),
321    problem_to_optimize_(NULL),
322    feasibility_mode_(false),
323    app_(NULL),
324    debug_apps_(),
325    testOthers_(false),
326    warmstart_(NULL),
327    rowsense_(NULL),
328    rhs_(NULL),
329    rowrange_(NULL),
330    reducedCosts_(NULL),
331    OsiDualObjectiveLimit_(1e200),
332    hasVarNamesFile_(true),
333    nCallOptimizeTNLP_(0),
334    totalNlpSolveTime_(0),
335    totalIterations_(0),
336    maxRandomRadius_(1e08),
337    randomGenerationType_(0),
338    max_perturbation_(COIN_DBL_MAX),
339    pushValue_(1e-02),
340    numRetryInitial_(-1),
341    numRetryResolve_(-1),
342    numRetryInfeasibles_(-1),
343    numRetryUnsolved_(1),
344    messages_(),
345    pretendFailIsInfeasible_(0),
346    hasContinuedAfterNlpFailure_(false),
347    numIterationSuspect_(-1),
348    hasBeenOptimized_(false),
349    obj_(NULL),
350    feasibilityProblem_(NULL),
351    jRow_(NULL),
352    jCol_(NULL),
353    jValues_(NULL),
354    nnz_jac(0),
355    constTypes_(NULL),
356    nNonLinear_(0),
357    tiny_(1e-8),
358    veryTiny_(1e-20),
359    infty_(1e100),
360    exposeWarmStart_(false),
361    firstSolve_(true),
362    cutStrengthener_(NULL),
363    oaMessages_(),
364    oaHandler_(NULL)
365{
366   oaHandler_ = new OaMessageHandler;
367   oaHandler_->setLogLevel(0);
368}
369
370void 
371OsiTMINLPInterface::initialize(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
372                               Ipopt::SmartPtr<Ipopt::OptionsList> options,
373                               Ipopt::SmartPtr<Ipopt::Journalist> journalist,
374                               const std::string & prefix,
375                               Ipopt::SmartPtr<TMINLP> tminlp){
376  if(!IsValid(app_))
377     createApplication(roptions, options, journalist, prefix);
378  setModel(tminlp); 
379}
380
381void OsiTMINLPInterface::setSolver(Ipopt::SmartPtr<TNLPSolver> app){
382  app_ = app;
383  }
384
385
386void
387OsiTMINLPInterface::createApplication(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
388                                      Ipopt::SmartPtr<Ipopt::OptionsList> options,
389                                      Ipopt::SmartPtr<Ipopt::Journalist> journalist,
390                                      const std::string & prefix
391                                      )
392{
393  assert(!IsValid(app_));
394  int ival;
395  options->GetEnumValue("nlp_solver", ival, prefix.c_str());
396  Solver s = (Solver) ival;
397  if(s == EFilterSQP){
398    testOthers_ = false;;
399#ifdef COIN_HAS_FILTERSQP
400    app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
401#else
402   throw SimpleError("createApplication",
403                     "Bonmin not configured to run with FilterSQP.");
404#endif   
405
406#ifdef COIN_HAS_IPOPT
407   debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix)); 
408#endif
409  }
410  else if(s == EIpopt){
411    testOthers_ = false;
412#ifdef COIN_HAS_IPOPT
413    app_ = new IpoptSolver(roptions, options, journalist, prefix);
414#else
415   throw SimpleError("createApplication",
416                     "Bonmin not configured to run with Ipopt.");
417#endif
418#ifdef COIN_HAS_FILTERSQP
419    debug_apps_.push_back(new Bonmin::FilterSolver(roptions, options, journalist, prefix));
420#endif
421  }
422  else if(s == EAll){
423#ifdef COIN_HAS_FILTERSQP
424    app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
425#else
426   throw SimpleError("createApplication",
427                     "Bonmin not configured to run with Ipopt.");
428#endif
429#ifdef COIN_HAS_IPOPT
430   debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix)); 
431#endif
432    testOthers_ = true;
433  }
434  if (!app_->Initialize("")) {
435    throw CoinError("Error during initialization of app_","createApplication", "OsiTMINLPInterface");
436  }
437  for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin() ; 
438      i != debug_apps_.end() ; i++){
439    (*i)->Initialize("");
440  }
441  extractInterfaceParams();
442 
443}
444
445/** Facilitator to allocate a tminlp and an application. */
446void
447OsiTMINLPInterface::setModel(SmartPtr<TMINLP> tminlp)
448{
449  assert(IsValid(tminlp));
450  tminlp_ = tminlp;
451  problem_ = new TMINLP2TNLP(tminlp_);
452  feasibilityProblem_ = new TNLP2FPNLP
453        (SmartPtr<TNLP>(GetRawPtr(problem_)));
454  if(feasibility_mode_){
455    problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
456  }
457  else {
458    problem_to_optimize_ = GetRawPtr(problem_);
459  }
460}
461
462
463
464void
465OsiTMINLPInterface::readOptionFile(const std::string & fileName)
466{
467  if(IsValid(app_)){
468  std::ifstream is;
469  if (fileName != "") {
470    try {
471      is.open(fileName.c_str());
472    }
473    catch(std::bad_alloc) {
474      throw CoinError("Not enough memory to open option file.\n", "readOptionFile", "OsiTMINLPInterface");
475    }
476  }
477  options()->ReadFromStream(*app_->journalist(), is);
478  extractInterfaceParams();
479  }
480}
481
482/// Copy constructor
483OsiTMINLPInterface::OsiTMINLPInterface (const OsiTMINLPInterface &source):
484    OsiSolverInterface(source),
485    tminlp_(source.tminlp_),
486    problem_(NULL),
487    problem_to_optimize_(NULL),
488    feasibility_mode_(source.feasibility_mode_),
489    rowsense_(NULL),
490    rhs_(NULL),
491    rowrange_(NULL),
492    reducedCosts_(NULL),
493    OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
494    hasVarNamesFile_(source.hasVarNamesFile_),
495    nCallOptimizeTNLP_(0),
496    totalNlpSolveTime_(0),
497    totalIterations_(0),
498    maxRandomRadius_(source.maxRandomRadius_),
499    randomGenerationType_(source.randomGenerationType_),
500    max_perturbation_(source.max_perturbation_),
501    pushValue_(source.pushValue_),
502    numRetryInitial_(source.numRetryInitial_),
503    numRetryResolve_(source.numRetryResolve_),
504    numRetryInfeasibles_(source.numRetryInfeasibles_),
505    numRetryUnsolved_(source.numRetryUnsolved_),
506    messages_(),
507    pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
508    hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
509    numIterationSuspect_(source.numIterationSuspect_),
510    hasBeenOptimized_(source.hasBeenOptimized_),
511    obj_(NULL),
512    feasibilityProblem_(NULL),
513    jRow_(NULL),
514    jCol_(NULL),
515    jValues_(NULL),
516    nnz_jac(source.nnz_jac),
517    constTypes_(NULL),
518//    constTypesNum_(NULL),
519    nNonLinear_(0),
520    tiny_(source.tiny_),
521    veryTiny_(source.veryTiny_),
522    infty_(source.infty_),
523    exposeWarmStart_(source.exposeWarmStart_),
524    firstSolve_(true),
525    cutStrengthener_(source.cutStrengthener_),
526    oaMessages_(),
527    oaHandler_(NULL),
528    strong_branching_solver_(source.strong_branching_solver_)
529{
530   if(0 && defaultHandler()){
531     messageHandler()->setLogLevel(source.messageHandler()->logLevel());
532   }
533  //Pass in message handler
534  if(0 && source.messageHandler())
535    passInMessageHandler(source.messageHandler());
536   //Copy options from old application
537  if(IsValid(source.tminlp_)) {
538    problem_ = source.problem_->clone();
539    feasibilityProblem_ = new TNLP2FPNLP
540        (SmartPtr<TNLP>(GetRawPtr(problem_)), source.feasibilityProblem_);
541    if(feasibility_mode_)
542      problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
543    else
544      problem_to_optimize_ = GetRawPtr(problem_);
545    pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
546
547    setAuxiliaryInfo(source.getAuxiliaryInfo());
548    // Copy options from old application
549    app_ = source.app_->clone();
550    for(std::list<Ipopt::SmartPtr<TNLPSolver> >::const_iterator i = source.debug_apps_.begin();
551        i != source.debug_apps_.end() ; i++){
552      debug_apps_.push_back((*i)->clone());
553    }
554    testOthers_ = source.testOthers_;
555  }
556  else {
557    throw SimpleError("Don't know how to copy an empty IpoptInterface.",
558        "copy constructor");
559  }
560
561  warmstart_ = source.warmstart_ ? source.warmstart_->clone() : NULL;
562
563  if(source.obj_) {
564    obj_ = new double[source.getNumCols()];
565    CoinCopyN(source.obj_, source.getNumCols(), obj_);
566  }
567
568
569   oaHandler_ = new OaMessageHandler(*source.oaHandler_);;
570
571}
572
573OsiSolverInterface * 
574OsiTMINLPInterface::clone(bool copyData ) const
575{
576  if(copyData)
577    return new OsiTMINLPInterface(*this);
578  else return new OsiTMINLPInterface;
579}
580
581/// Assignment operator
582OsiTMINLPInterface & OsiTMINLPInterface::operator=(const OsiTMINLPInterface& rhs)
583{
584  if(this!= &rhs) {
585    OsiSolverInterface::operator=(rhs);
586    OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
587    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
588    totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
589    totalIterations_ = rhs.totalIterations_;
590    maxRandomRadius_ = rhs.maxRandomRadius_;
591    hasVarNamesFile_ = rhs.hasVarNamesFile_;
592    pushValue_ = rhs.pushValue_;
593
594    delete warmstart_;
595    warmstart_ = NULL;
596
597    if(IsValid(rhs.tminlp_)) {
598
599      tminlp_ = rhs.tminlp_;
600      problem_ = new TMINLP2TNLP(tminlp_);
601      problem_to_optimize_ = GetRawPtr(problem_);
602      app_ = rhs.app_->clone();
603
604      warmstart_ = rhs.warmstart_ ? rhs.warmstart_->clone() : NULL;
605
606      feasibilityProblem_ = new TNLP2FPNLP
607          (SmartPtr<TNLP>(GetRawPtr(problem_)));
608      nnz_jac = rhs.nnz_jac;
609
610      if(constTypes_ != NULL) {
611        delete [] constTypes_;
612        constTypes_ = NULL;
613      }
614      if(rhs.constTypes_ != NULL) {
615        constTypes_ = new TNLP::LinearityType[getNumRows()];
616        CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
617      }
618/*
619      if(constTypesNum_ != NULL) {
620        delete [] constTypesNum_;
621        constTypesNum_ = NULL;
622      }
623      if(rhs.constTypesNum_ != NULL) {
624        constTypesNum_ = new int[getNumRows()];
625        CoinCopyN(rhs.constTypesNum_, getNumRows(), constTypesNum_);
626      }
627*/
628      if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
629        jValues_ = new double [nnz_jac];
630        jCol_    = new Index [nnz_jac];
631        jRow_    = new Index [nnz_jac];
632        CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
633        CoinCopyN(rhs.jCol_    , nnz_jac,jCol_    );
634        CoinCopyN(rhs.jRow_    , nnz_jac,jRow_    );
635      }
636      else if(nnz_jac > 0) {
637        throw CoinError("Arrays for storing jacobian are inconsistant.",
638            "copy constructor",
639            "IpoptOAInterface");
640      }
641      tiny_ = rhs.tiny_;
642      veryTiny_ = rhs.veryTiny_;
643      infty_ = rhs.infty_;
644      exposeWarmStart_ = rhs.exposeWarmStart_;
645
646    }
647    else {
648      tminlp_ =NULL;
649      problem_ = NULL;
650      app_ = NULL;
651      feasibilityProblem_ = NULL;
652    }
653
654
655    if(obj_) {
656      delete [] obj_;
657      obj_ = NULL;
658    }
659    if(rhs.obj_) {
660      obj_ = new double[rhs.getNumCols()];
661      CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
662    }
663
664    hasVarNamesFile_ = rhs.hasVarNamesFile_;
665
666    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
667    totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
668    totalIterations_ = rhs.totalIterations_;
669    maxRandomRadius_ = rhs.maxRandomRadius_;
670    pushValue_ = rhs.pushValue_;
671    numRetryInitial_ = rhs.numRetryInitial_;
672    numRetryResolve_ = rhs.numRetryResolve_;
673    numRetryInfeasibles_ = rhs.numRetryInfeasibles_;
674    numRetryUnsolved_ = rhs.numRetryUnsolved_;
675    pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
676    numIterationSuspect_ = rhs.numIterationSuspect_;
677
678    hasBeenOptimized_ = rhs.hasBeenOptimized_;
679    cutStrengthener_ = rhs.cutStrengthener_;
680
681    delete oaHandler_;
682    oaHandler_ = new OaMessageHandler(*rhs.oaHandler_);
683    strong_branching_solver_ = rhs.strong_branching_solver_;
684
685    freeCachedData();
686  }
687  return *this;
688}
689
690const SmartPtr<OptionsList> OsiTMINLPInterface::options() const
691{
692  if(!IsValid(app_)) {
693    messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
694    return NULL;
695  }
696  else
697    return app_->options();
698}
699
700SmartPtr<OptionsList> OsiTMINLPInterface::options()
701{
702  if(!IsValid(app_)) {
703    messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
704    return NULL;
705  }
706  else
707    return app_->options();
708}
709
710/// Destructor
711OsiTMINLPInterface::~OsiTMINLPInterface ()
712{
713  freeCachedData();
714  delete [] jRow_;
715  delete [] jCol_;
716  delete [] jValues_;
717  delete [] constTypes_;
718  delete [] obj_;
719  delete oaHandler_;
720  delete warmstart_;
721}
722
723void
724OsiTMINLPInterface::freeCachedColRim()
725{
726    if(reducedCosts_!=NULL) {
727    delete []  reducedCosts_;
728    reducedCosts_ = NULL;
729  }
730
731}
732
733void
734OsiTMINLPInterface::freeCachedRowRim()
735{
736  if(rowsense_!=NULL) {
737    delete []  rowsense_;
738    rowsense_ = NULL;
739  }
740  if(rhs_!=NULL) {
741    delete []  rhs_;
742    rhs_ = NULL;
743  }
744  if(rowrange_!=NULL) {
745    delete []  rowrange_;
746    rowrange_ = NULL;
747  }
748  //   if(dualsol_!=NULL)
749  //     {
750  //       delete []  dualsol_; dualsol_ = NULL;
751  //     }
752}
753
754void
755OsiTMINLPInterface::freeCachedData()
756{
757  freeCachedColRim();
758  freeCachedRowRim();
759}
760
761const char * OsiTMINLPInterface::OPT_SYMB="OPT";
762const char * OsiTMINLPInterface::FAILED_SYMB="FAILED";
763const char * OsiTMINLPInterface::UNBOUND_SYMB="UNBOUNDED";
764const char * OsiTMINLPInterface::INFEAS_SYMB="INFEAS";
765const char * OsiTMINLPInterface::TIME_SYMB="TIME";
766///////////////////////////////////////////////////////////////////
767// WarmStart Information                                                                           //
768///////////////////////////////////////////////////////////////////
769
770
771void
772OsiTMINLPInterface::resolveForCost(int numsolve, bool keepWarmStart)
773{
774  // This method assumes that a problem has just been solved and we try for a
775  // different solution. So disregard (in fact, clear out) any warmstart info
776  // we might have, and acquire a new one before returning.
777  delete warmstart_;
778  warmstart_ = NULL;
779 
780  Coin::SmartPtr<SimpleReferencedPtr<CoinWarmStart> > ws_backup = NULL;
781  if(!exposeWarmStart_ && keepWarmStart){
782    //if warm start is not exposed, we need to store the current starting point to
783    //restore it at the end of the method
784    ws_backup = make_referenced(app_->getUsedWarmStart(problem_));
785  }
786  /** Save current optimum. */
787  vector<double> point(getNumCols()*3+ getNumRows());
788  double bestBound = getObjValue();
789  CoinCopyN(getColSolution(),
790      getNumCols(), point());
791  CoinCopyN(getRowPrice(),
792      2*getNumCols()+ getNumRows(),
793      point() + getNumCols());
794
795  if(isProvenOptimal())
796    messageHandler()->message(SOLUTION_FOUND,
797        messages_)
798    <<1<<getObjValue()<<bestBound
799    <<CoinMessageEol;
800  else
801    messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
802        messages_)
803    <<1
804    <<CoinMessageEol;
805  for(int f = 0; f < numsolve ; f++) {
806    messageHandler()->message(WARNING_RESOLVING,
807        messages_)
808    <<f+1<< CoinMessageEol ;
809    randomStartingPoint();
810    solveAndCheckErrors(0,0,"resolve cost");
811
812
813    char c=' ';
814    //Is solution better than previous
815    if(isProvenOptimal() &&
816        getObjValue()<bestBound) {
817      c='*';
818      messageHandler()->message(BETTER_SOL, messages_)<<getObjValue()<<f+1<< CoinMessageEol;
819      CoinCopyN(getColSolution(),
820          getNumCols(), point());
821      CoinCopyN(getRowPrice(),
822          2*getNumCols()+ getNumRows(),
823          point() + getNumCols());
824      bestBound = getObjValue();
825    }
826
827    messageHandler()->message(LOG_LINE, messages_)
828    <<c<<f+1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<app_->CPUTime()<<"resolve cost"<<CoinMessageEol;
829
830
831    if(isProvenOptimal())
832      messageHandler()->message(SOLUTION_FOUND,
833          messages_)
834      <<f+2<<getObjValue()<<bestBound
835      <<CoinMessageEol;
836    else if(!isAbandoned())
837      messageHandler()->message(UNSOLVED_PROBLEM_FOUND,
838          messages_)
839      <<f+2
840      <<CoinMessageEol;
841    else
842      messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
843          messages_)
844      <<f+2
845      <<CoinMessageEol;
846  }
847  setColSolution(point());
848  setRowPrice(point() + getNumCols());
849  app_->enableWarmStart();
850
851  optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
852  hasBeenOptimized_ = true;
853
854  if(!exposeWarmStart_ && keepWarmStart) {
855    app_->setWarmStart(ws_backup->ptr(), problem_);
856  }
857}
858
859void
860OsiTMINLPInterface::resolveForRobustness(int numsolve)
861{
862  // This method assumes that a problem has just been solved and we try for a
863  // different solution. So disregard (in fact, clear out) any warmstart info
864  // we might have, and acquire a new one before returning.
865  delete warmstart_;
866  warmstart_ = NULL;
867 
868
869  CoinWarmStart * ws_backup = NULL;
870  if(!exposeWarmStart_){
871    //if warm start is not exposed, we need to store the current starting point to
872    //restore it at the end of the method
873    ws_backup = app_->getUsedWarmStart(problem_);
874  }
875  //std::cerr<<"Resolving the problem for robustness"<<std::endl;
876  //First remove warm start point and resolve
877  app_->disableWarmStart();
878  problem()->resetStartingPoint();
879  messageHandler()->message(WARNING_RESOLVING,
880      messages_)
881  <<1<< CoinMessageEol ;
882  solveAndCheckErrors(0,0,"resolve robustness");
883
884
885  char c='*';
886  if(isAbandoned()) {
887    c=' ';
888  }
889  messageHandler()->message(LOG_LINE, messages_)
890  <<c<<1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<
891    app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
892
893
894  if(!isAbandoned()) {
895    messageHandler()->message(WARN_SUCCESS_WS, messages_) << CoinMessageEol ;
896    // re-enable warmstart and get it
897    app_->enableWarmStart();
898    if (!exposeWarmStart_) {
899      app_->setWarmStart(ws_backup, problem_);
900      delete ws_backup;
901    }
902    return; //we won go on
903  }
904
905  //still unsolved try again with different random starting points
906  for(int f = 0; f < numsolve ; f++) {
907    messageHandler()->message(WARNING_RESOLVING,
908        messages_)
909    <<f+2<< CoinMessageEol ;
910
911    randomStartingPoint();
912    solveAndCheckErrors(0,0,"resolve robustness");
913
914
915    messageHandler()->message(IPOPT_SUMMARY, messages_)
916    <<"resolveForRobustness"<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
917
918
919    char c='*';
920    if(isAbandoned()) {
921      c=' ';
922    }
923    messageHandler()->message(LOG_LINE, messages_)
924    <<c<<f+2<<statusAsString()<<getObjValue()
925    <<app_->IterationCount()<<app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
926
927
928    if(!isAbandoned()) {
929      messageHandler()->message(WARN_SUCCESS_RANDOM, messages_)
930        << f+2 << CoinMessageEol ;
931      // re-enable warmstart and get it
932      app_->enableWarmStart();
933      if (!exposeWarmStart_) {
934        app_->setWarmStart(ws_backup, problem_);
935        delete ws_backup;
936      }
937       
938      return; //we have found a solution and continue
939    }
940  }
941
942
943  if(!exposeWarmStart_){
944    app_->setWarmStart(ws_backup, problem_);
945    delete ws_backup;
946  }
947  if(pretendFailIsInfeasible_) {
948    if(pretendFailIsInfeasible_ == 1) {
949      messageHandler()->message(WARN_CONTINUING_ON_FAILURE,
950          messages_)
951      <<CoinMessageEol;
952      hasContinuedAfterNlpFailure_ = 1;
953    }
954    return;
955  }
956  else {
957    std::string probName;
958    getStrParam(OsiProbName,probName);
959    throw newUnsolvedError(app_->errorCode(), problem_,
960                           probName);
961  }
962  // Do NOT get warmstart in other cases
963}
964
965////////////////////////////////////////////////////////////////////
966// Problem information methods                                    //
967////////////////////////////////////////////////////////////////////
968/// Get number of columns
969int OsiTMINLPInterface::getNumCols() const
970{
971
972  return problem_->num_variables();
973}
974
975
976/// Get number of rows
977int
978OsiTMINLPInterface::getNumRows() const
979{
980  return problem_->num_constraints();
981}
982
983const double *
984OsiTMINLPInterface::getColLower() const
985{
986  return problem_->x_l();
987}
988
989const double *
990OsiTMINLPInterface::getColUpper() const
991{
992  return problem_->x_u();
993}
994
995#if 1
996
997
998///get name of variables
999const OsiSolverInterface::OsiNameVec& 
1000OsiTMINLPInterface::getVarNames() {
1001  return getColNames();
1002}
1003#endif
1004
1005
1006void OsiTMINLPInterface::extractSenseRhsAndRange() const
1007{
1008  assert(rowsense_==NULL&&rhs_==NULL&&rowrange_==NULL);
1009  int numrows = problem_->num_constraints();
1010  if(numrows == 0) return;
1011  const double * rowLower = getRowLower();
1012  const double * rowUpper = getRowUpper();
1013  rowsense_ = new char [numrows];
1014  rhs_ = new double [numrows];
1015  rowrange_ = new double [numrows];
1016  for(int i = 0 ; i < numrows ; i++) {
1017    rowrange_[i]=0.;
1018    convertBoundToSense(rowLower[i], rowUpper[i], rowsense_[i], rhs_[i], rowrange_[i]);
1019  }
1020}
1021
1022/** Get pointer to array[getNumRows()] of row constraint senses.
1023    <ul>
1024    <li>'L': <= constraint
1025    <li>'E': =  constraint
1026    <li>'G': >= constraint
1027    <li>'R': ranged constraint
1028    <li>'N': free constraint
1029    </ul>
1030*/
1031const char *
1032OsiTMINLPInterface::getRowSense() const
1033{
1034  if(rowsense_==NULL) {
1035    extractSenseRhsAndRange();
1036  }
1037  return rowsense_;
1038}
1039
1040/** Get pointer to array[getNumRows()] of rows right-hand sides
1041    <ul>
1042    <li> if rowsense()[i] == 'L' then rhs()[i] == rowupper()[i]
1043    <li> if rowsense()[i] == 'G' then rhs()[i] == rowlower()[i]
1044    <li> if rowsense()[i] == 'R' then rhs()[i] == rowupper()[i]
1045    <li> if rowsense()[i] == 'N' then rhs()[i] == 0.0
1046    </ul>
1047*/
1048const double *
1049OsiTMINLPInterface::getRightHandSide() const
1050{
1051  if(rhs_==NULL) {
1052    extractSenseRhsAndRange();
1053  }
1054  return rhs_;
1055}
1056
1057/** Get pointer to array[getNumRows()] of row ranges.
1058    <ul>
1059    <li> if rowsense()[i] == 'R' then
1060    rowrange()[i] == rowupper()[i] - rowlower()[i]
1061    <li> if rowsense()[i] != 'R' then
1062    rowrange()[i] is 0.0
1063    </ul>
1064*/
1065const double *
1066OsiTMINLPInterface::getRowRange() const
1067{
1068  if(rowrange_==NULL) {
1069    extractSenseRhsAndRange();
1070  }
1071  return rowrange_;
1072}
1073
1074const double *
1075OsiTMINLPInterface::getRowLower() const
1076{
1077  return problem_->g_l();
1078}
1079
1080const double *
1081OsiTMINLPInterface::getRowUpper() const
1082{
1083  return problem_->g_u();
1084}
1085
1086/// Return true if column is continuous
1087bool
1088OsiTMINLPInterface::isContinuous(int colNumber) const
1089{
1090  return (problem_->var_types()[colNumber]==TMINLP::CONTINUOUS);
1091}
1092
1093/// Return true if column is binary
1094bool
1095OsiTMINLPInterface::isBinary(int colNumber) const
1096{
1097  return (problem_->var_types()[colNumber]==TMINLP::BINARY);
1098}
1099
1100/** Return true if column is integer.
1101    Note: This function returns true if the the column
1102    is binary or a general integer.
1103*/
1104bool
1105OsiTMINLPInterface::isInteger(int colNumber) const
1106{
1107  return ((problem_->var_types()[colNumber]==TMINLP::BINARY)||
1108      (problem_->var_types()[colNumber]==TMINLP::INTEGER));
1109}
1110
1111/// Return true if column is general integer
1112bool
1113OsiTMINLPInterface::isIntegerNonBinary(int colNumber) const
1114{
1115  return (problem_->var_types()[colNumber]==TMINLP::INTEGER);
1116}
1117/// Return true if column is binary and not fixed at either bound
1118bool
1119OsiTMINLPInterface::isFreeBinary(int colNumber) const
1120{
1121  return ((problem_->var_types()[colNumber]==TMINLP::BINARY)
1122      &&((getColUpper()[colNumber]-getColLower()[colNumber]) > 1 - 1e-09));
1123}
1124
1125/// Get solver's value for infinity
1126double
1127OsiTMINLPInterface::getInfinity() const
1128{
1129  return COIN_DBL_MAX;
1130}
1131
1132/// Get pointer to array[getNumCols()] of primal solution vector
1133const double *
1134OsiTMINLPInterface::getColSolution() const
1135{
1136  if(hasBeenOptimized_)
1137    return problem_->x_sol();
1138  else
1139    return problem_->x_init();
1140}
1141
1142/// Get pointer to array[getNumRows()] of dual prices
1143const double *
1144OsiTMINLPInterface::getRowPrice() const
1145{
1146  if(hasBeenOptimized_)
1147    return problem_->duals_sol();
1148  else
1149    return problem_->duals_init();
1150}
1151
1152/// Get a pointer to array[getNumCols()] of reduced costs
1153const double *
1154OsiTMINLPInterface::getReducedCost() const
1155{
1156  (*handler_)<<"WARNING : trying to access reduced cost in Ipopt always retrun 0"<<CoinMessageEol;
1157  if(reducedCosts_==NULL) {
1158    reducedCosts_ = new double [getNumCols()];
1159    CoinFillN(reducedCosts_,getNumCols(),0.);
1160  }
1161  return reducedCosts_;
1162}
1163
1164/** Get pointer to array[getNumRows()] of row activity levels (constraint
1165    matrix times the solution vector */
1166const double *
1167OsiTMINLPInterface::getRowActivity() const
1168{
1169  return problem_->g_sol();
1170}
1171
1172/** Get how many iterations it took to solve the problem (whatever
1173    "iteration" mean to the solver.
1174*/
1175int
1176OsiTMINLPInterface::getIterationCount() const
1177{
1178    return app_->IterationCount();
1179}
1180
1181
1182/** Set a single column lower bound.
1183    Use -getInfinity() for -infinity. */
1184void
1185OsiTMINLPInterface::setColLower( int elementIndex, double elementValue )
1186{
1187  //  if(fabs(problem_->x_l()[elementIndex]-elementValue)>1e-06)
1188  problem_->SetVariableLowerBound(elementIndex,elementValue);
1189  hasBeenOptimized_ = false;
1190}
1191
1192/** Set a single column upper bound.
1193    Use getInfinity() for infinity. */
1194void
1195OsiTMINLPInterface::setColUpper( int elementIndex, double elementValue )
1196{
1197  //  if(fabs(problem_->x_u()[elementIndex]-elementValue)>1e-06)
1198  problem_->SetVariableUpperBound(elementIndex,elementValue);
1199  hasBeenOptimized_ = false;
1200}
1201
1202/** Set the lower bounds for all columns
1203    Use -getInfinity() for -infinity. */
1204void
1205OsiTMINLPInterface::setColLower( const double* array )
1206{
1207  problem_->SetVariablesLowerBounds(problem_->num_variables(),
1208                                  array);
1209  hasBeenOptimized_ = false;
1210}
1211
1212/** Set Set the upper bounds for all columns
1213    Use getInfinity() for infinity. */
1214void
1215OsiTMINLPInterface::setColUpper( const double* array )
1216{
1217  problem_->SetVariablesUpperBounds(problem_->num_variables(), 
1218                                  array);
1219  hasBeenOptimized_ = false;
1220}
1221
1222/** Set a single row lower bound.
1223    Use -getInfinity() for -infinity. */
1224void
1225OsiTMINLPInterface::setRowLower( int elementIndex, double elementValue )
1226{
1227  throw SimpleError("Not implemented yet but should be if necessary.",
1228      "setRowLower");
1229  hasBeenOptimized_ = false;
1230}
1231
1232/** Set a single row upper bound.
1233    Use getInfinity() for infinity. */
1234void
1235OsiTMINLPInterface::setRowUpper( int elementIndex, double elementValue )
1236{
1237  throw SimpleError("Not implemented yet but should be if necessary.",
1238      "setRowUpper");
1239  hasBeenOptimized_ = false;
1240}
1241
1242/** Set the type of a single row */
1243void
1244OsiTMINLPInterface::setRowType(int index, char sense, double rightHandSide,
1245    double range)
1246{
1247  throw SimpleError("Not implemented yet but should be if necessary.",
1248      "setRowType");
1249  hasBeenOptimized_ = false;
1250}
1251
1252
1253/// Set the objective function sense.
1254/// (1 for min (default), -1 for max)
1255void
1256OsiTMINLPInterface::setObjSense(double s)
1257{
1258  throw SimpleError("Can not change objective sense of an Ipopt problem.",
1259      "setObjSense");
1260  hasBeenOptimized_ = false;
1261}
1262
1263/** Set the primal solution variable values
1264
1265colsol[getNumCols()] is an array of values for the primal variables.
1266These values are copied to memory owned by the solver interface object
1267or the solver.  They will be returned as the result of getColSolution()
1268until changed by another call to setColSolution() or by a call to any
1269solver routine.  Whether the solver makes use of the solution in any
1270way is solver-dependent.
1271*/
1272void
1273OsiTMINLPInterface::setColSolution(const double *colsol)
1274{
1275  problem_->setxInit(getNumCols(), colsol);
1276  hasBeenOptimized_ = false;
1277}
1278
1279/** Set dual solution variable values
1280
1281rowprice[getNumRows()] is an array of values for the dual
1282variables. These values are copied to memory owned by the solver
1283interface object or the solver.  They will be returned as the result of
1284getRowPrice() until changed by another call to setRowPrice() or by a
1285call to any solver routine.  Whether the solver makes use of the
1286solution in any way is solver-dependent.
1287*/
1288
1289void
1290OsiTMINLPInterface::setRowPrice(const double * rowprice)
1291{
1292  problem_->setDualsInit(getNumCols()*2 + getNumRows(), rowprice);
1293  hasBeenOptimized_ = false;
1294}
1295
1296  /*! \brief Get an empty warm start object
1297
1298  This routine returns an empty CoinWarmStartBasis object. Its purpose is
1299  to provide a way to give a client a warm start basis object of the
1300  appropriate type, which can resized and modified as desired.
1301  */
1302CoinWarmStart *
1303OsiTMINLPInterface::getEmptyWarmStart () const
1304{return app_->getEmptyWarmStart();}
1305
1306  /** Get warmstarting information */
1307CoinWarmStart* 
1308OsiTMINLPInterface::getWarmStart() const
1309{
1310  if (exposeWarmStart_) {
1311    return internal_getWarmStart();;
1312  }
1313  else {
1314    return getEmptyWarmStart();
1315  }
1316}
1317  /** Set warmstarting information. Return true/false depending on whether
1318      the warmstart information was accepted or not. */
1319bool 
1320OsiTMINLPInterface::setWarmStart(const CoinWarmStart* ws)
1321{
1322  if (exposeWarmStart_) {
1323    return internal_setWarmStart(ws);
1324  }
1325  else {
1326    return true;
1327  }
1328}
1329  /** Get warmstarting information */
1330CoinWarmStart* 
1331OsiTMINLPInterface::internal_getWarmStart() const
1332{
1333  if (exposeWarmStart_ && warmstart_) {
1334    return warmstart_->clone();
1335  }
1336  else {
1337    return getEmptyWarmStart();
1338  }
1339}
1340  /** Set warmstarting information. Return true/false depending on whether
1341      the warmstart information was accepted or not. */
1342bool 
1343OsiTMINLPInterface::internal_setWarmStart(const CoinWarmStart* ws)
1344{
1345  delete warmstart_;
1346  warmstart_ = NULL;
1347  hasBeenOptimized_ = false;
1348  if (exposeWarmStart_) {
1349    if (ws == NULL) {
1350      return true;
1351    }
1352  if(app_->warmStartIsValid(ws)) {
1353    warmstart_ = ws->clone();
1354    return true;
1355  }
1356  // See if it is anything else than the CoinWarmStartBasis that all others
1357  // derive from
1358  const CoinWarmStartPrimalDual* pdws =
1359    dynamic_cast<const CoinWarmStartPrimalDual*>(ws);
1360  if (pdws) {
1361    // Must be an IpoptWarmStart, since the others do not derive from this.
1362    // Accept it.
1363    warmstart_ = new IpoptWarmStart(*pdws);
1364    return true;
1365  }
1366  return false;
1367  }
1368  else {
1369    return true;
1370  }
1371}
1372
1373/** Set the index-th variable to be a continuous variable */
1374void
1375OsiTMINLPInterface::setContinuous(int index)
1376{
1377  problem_->SetVariableType(index, TMINLP::CONTINUOUS);
1378  hasBeenOptimized_ = false;
1379}
1380/** Set the index-th variable to be an integer variable */
1381void
1382OsiTMINLPInterface::setInteger(int index)
1383{
1384  problem_->SetVariableType(index, TMINLP::INTEGER);
1385  hasBeenOptimized_ = false;
1386}
1387
1388/// Get objective function value (can't use default)
1389double
1390OsiTMINLPInterface::getObjValue() const
1391{
1392  return problem_->obj_value();
1393}
1394
1395//#############################################################################
1396// Parameter related methods
1397//#############################################################################
1398
1399bool
1400OsiTMINLPInterface::setIntParam(OsiIntParam key, int value)
1401{
1402  //  debugMessage("OsiCpxSolverInterface::setIntParam(%d, %d)\n", key, value);
1403
1404  bool retval = false;
1405  switch (key) {
1406  case OsiMaxNumIteration:
1407    retval = false;
1408    break;
1409  case OsiMaxNumIterationHotStart:
1410    if( value >= 0 ) {
1411      retval = false;
1412    }
1413    else
1414      retval = false;
1415    break;
1416  case OsiLastIntParam:
1417    retval = false;
1418    break;
1419  default:
1420    retval = false;
1421    (*handler_)<< "Unhandled case in setIntParam\n"<<CoinMessageEol;
1422    break;
1423  }
1424  return retval;
1425}
1426
1427//-----------------------------------------------------------------------------
1428
1429bool
1430OsiTMINLPInterface::setDblParam(OsiDblParam key, double value)
1431{
1432  //  debugMessage("OsiTMINLPInterface::setDblParam(%d, %g)\n", key, value);
1433
1434  bool retval = false;
1435  switch (key) {
1436  case OsiDualObjectiveLimit:
1437    OsiDualObjectiveLimit_ = value;
1438    retval = true;
1439    break;
1440  case OsiPrimalObjectiveLimit:
1441    (*handler_)<<"Can not set primal objective limit parameter"<<CoinMessageEol;
1442    retval = false;
1443    break;
1444  case OsiDualTolerance:
1445    (*handler_)<<"Can not set dual tolerance parameter"<<CoinMessageEol;
1446    retval = false;
1447    break;
1448  case OsiPrimalTolerance:
1449    (*handler_)<<"Can not set primal tolerance parameter"<<CoinMessageEol;
1450    retval = false;
1451  case OsiObjOffset:
1452    retval = OsiSolverInterface::setDblParam(key,value);
1453    break;
1454  case OsiLastDblParam:
1455    retval = false;
1456    break;
1457  default:
1458    retval = false;
1459    (*handler_) << "Unhandled case in setDblParam"<<CoinMessageEol;
1460    break;
1461  }
1462  return retval;
1463}
1464
1465
1466//-----------------------------------------------------------------------------
1467
1468bool
1469OsiTMINLPInterface::setStrParam(OsiStrParam key, const std::string & value)
1470{
1471  //  debugMessage("OsiTMINLPInterface::setStrParam(%d, %s)\n", key, value.c_str());
1472
1473  bool retval=false;
1474  switch (key) {
1475  case OsiProbName:
1476    OsiSolverInterface::setStrParam(key,value);
1477    return retval = true;
1478  case OsiSolverName:
1479    return false;
1480  case OsiLastStrParam:
1481    return false;
1482  }
1483  return false;
1484}
1485
1486//-----------------------------------------------------------------------------
1487
1488bool
1489OsiTMINLPInterface::getIntParam(OsiIntParam key, int& value) const
1490{
1491  //  debugMessage("OsiTMINLPInterface::getIntParam(%d)\n", key);
1492
1493  value = -COIN_INT_MAX; // Give a dummy value
1494  bool retval = false;
1495  switch (key) {
1496  case OsiMaxNumIteration:
1497    retval = false;
1498    break;
1499  case OsiMaxNumIterationHotStart:
1500    retval = false;
1501    break;
1502  case OsiLastIntParam:
1503    retval = false;
1504    break;
1505  default:
1506    retval = false;
1507    (*handler_) << "Unhandled case in setIntParam"<<CoinMessageEol;
1508  }
1509  return retval;
1510}
1511
1512//-----------------------------------------------------------------------------
1513
1514bool
1515OsiTMINLPInterface::getDblParam(OsiDblParam key, double& value) const
1516{
1517  //  debugMessage("OsiTMINLPInterface::getDblParam(%d)\n", key);
1518
1519  bool retval = false;
1520  switch (key) {
1521  case OsiDualObjectiveLimit:
1522    value = OsiDualObjectiveLimit_;
1523    retval = true;
1524    break;
1525  case OsiPrimalObjectiveLimit:
1526    value = getInfinity();
1527    retval = true;
1528    break;
1529  case OsiDualTolerance:
1530    retval = false;
1531    break;
1532  case OsiPrimalTolerance:
1533    options()->GetNumericValue("tol", value,"");
1534    value = 1e-07;
1535    retval = true;
1536    break;
1537  case OsiObjOffset:
1538    retval = OsiSolverInterface::getDblParam(key, value);
1539    break;
1540  case OsiLastDblParam:
1541    retval = false;
1542    break;
1543  }
1544  return retval;
1545}
1546
1547
1548//-----------------------------------------------------------------------------
1549
1550bool
1551OsiTMINLPInterface::getStrParam(OsiStrParam key, std::string & value) const
1552{
1553  //  debugMessage("OsiTMINLPInterface::getStrParam(%d)\n", key);
1554
1555  switch (key) {
1556  case OsiProbName:
1557    OsiSolverInterface::getStrParam(key, value);
1558    break;
1559  case OsiSolverName:
1560    value = "Ipopt";
1561    break;
1562  case OsiLastStrParam:
1563    return false;
1564  }
1565
1566  return true;
1567}
1568
1569void
1570OsiTMINLPInterface::randomStartingPoint()
1571{
1572  int numcols = getNumCols();
1573  const double * colLower = getColLower();
1574  const double * colUpper = getColUpper();
1575  double * sol = new double[numcols];
1576  const Number * x_init = problem_->x_init_user();
1577  const double* perturb_radius = NULL;
1578  if (randomGenerationType_ == perturb_suffix) {
1579    const TMINLP::PerturbInfo* pertubinfo = tminlp_->perturbInfo();
1580    if (pertubinfo) {
1581      perturb_radius = pertubinfo->GetPerturbationArray();
1582    }
1583    if (!perturb_radius) {
1584      throw SimpleError("Can't use perturb_radius if no radii are given.",
1585                        "randomStartingPoint");
1586    }
1587  }
1588  for(int i = 0 ; i < numcols ; i++) {
1589    int randomGenerationType = randomGenerationType_;
1590    if(x_init[i] < colLower[i] || x_init[i] > colUpper[i])
1591      randomGenerationType = uniform;
1592    if(randomGenerationType_ == uniform){
1593      double lower = std::min(-maxRandomRadius_,colUpper[i] - maxRandomRadius_);
1594      lower = std::max(colLower[i], lower);
1595      double upper = std::max(maxRandomRadius_,colLower[i] + maxRandomRadius_);
1596      upper = std::min(colUpper[i],upper);
1597      lower = std::min(upper,lower);
1598      upper = std::max(upper, lower);
1599      double interval = upper - lower;
1600      sol[i] = CoinDrand48()*(interval) + lower;}
1601    else if (randomGenerationType_ == perturb){
1602      const double lower = std::max(x_init[i] - max_perturbation_, colLower[i]);
1603      const double upper = std::min(x_init[i] + max_perturbation_, colUpper[i]);
1604      const double interval = upper - lower;
1605      sol[i]  = lower + CoinDrand48()*(interval);
1606    }
1607    else if (randomGenerationType_ == perturb_suffix){
1608      const double radius = perturb_radius[i];
1609      const double lower = std::max(x_init[i] - radius*max_perturbation_, colLower[i]);
1610      const double upper = std::min(x_init[i] + radius*max_perturbation_, colUpper[i]);
1611      const double interval = upper - lower;
1612      sol[i]  = lower + CoinDrand48()*(interval);
1613    }
1614  }
1615  app_->disableWarmStart();
1616  setColSolution(sol);
1617  delete [] sol;
1618}
1619
1620
1621
1622/** This methods initialiaze arrays for storing the jacobian */
1623int OsiTMINLPInterface::initializeJacobianArrays()
1624{
1625  Index n, m, nnz_h_lag;
1626  TNLP::IndexStyleEnum index_style;
1627  problem_to_optimize_->get_nlp_info( n, m, nnz_jac, nnz_h_lag, index_style);
1628
1629  if(jRow_ != NULL) delete [] jRow_;
1630  if(jCol_ != NULL) delete [] jCol_;
1631  if(jValues_ != NULL) delete [] jValues_;
1632
1633  jRow_ = new Index[nnz_jac];
1634  jCol_ = new Index[nnz_jac];
1635  jValues_ = new Number[nnz_jac];
1636  problem_to_optimize_->eval_jac_g(n, NULL, 0, m, nnz_jac, jRow_, jCol_, NULL);
1637  if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
1638  {
1639    for(int i = 0 ; i < nnz_jac ; i++){
1640      jRow_[i]--;
1641      jCol_[i]--;
1642    }
1643  }
1644
1645  if(constTypes_ != NULL) delete [] constTypes_;
1646//  if(constTypesNum_ != NULL) delete [] constTypesNum_;
1647
1648  constTypes_ = new TNLP::LinearityType[getNumRows()];
1649  problem_to_optimize_->get_constraints_linearity(getNumRows(), constTypes_);
1650//  constTypesNum_ = new int[getNumRows()];
1651  for(int i = 0; i < getNumRows() ; i++) {
1652    if(constTypes_[i]==TNLP::NON_LINEAR) {
1653      //constTypesNum_[i] =
1654      nNonLinear_++;
1655    }
1656  }
1657  return nnz_jac;
1658}
1659
1660
1661double 
1662OsiTMINLPInterface::getConstraintsViolation(const double *x, double &obj)
1663{
1664  int numcols = getNumCols();
1665  int numrows = getNumRows();
1666  double * g = new double[numrows];
1667  tminlp_->eval_g(numcols, x, 1, numrows, g);
1668  const double * rowLower = getRowLower();
1669  const double * rowUpper = getRowUpper();
1670
1671
1672  double norm = 0;
1673  for(int i = 0; i< numrows ; i++) {
1674    if(!constTypes_ || constTypes_[i] == TNLP::NON_LINEAR) {
1675      double rowViolation = 0;
1676      if(rowLower[i] > -1e10)
1677         rowViolation = std::max(0.,rowLower[i] - g[i]);
1678
1679      if(rowUpper[i] < 1e10);
1680        rowViolation = std::max(rowViolation, g[i] - rowUpper[i]);
1681
1682      norm = rowViolation > norm ? rowViolation : norm;
1683    }
1684  }
1685  tminlp_->eval_f(numcols, x, 1, obj);
1686  delete [] g;
1687  return norm;
1688}
1689
1690/** get infinity norm of constraint violation of a point and objective error*/
1691double
1692OsiTMINLPInterface::getNonLinearitiesViolation(const double *x, const double obj)
1693{
1694  double f;
1695  double norm = getConstraintsViolation(x, f);
1696  assert((f - obj) > -1e-08);
1697  norm =  (f - obj) > norm ? f - obj : norm;
1698  return norm;
1699}
1700
1701
1702
1703//A procedure to try to remove small coefficients in OA cuts (or make it non small
1704static inline
1705bool cleanNnz(double &value, double colLower, double colUpper,
1706    double rowLower, double rowUpper, double colsol,
1707    double & lb, double &ub, double tiny, double veryTiny)
1708{
1709  if(fabs(value)>= tiny) return 1;
1710
1711  if(fabs(value)<veryTiny) return 0;//Take the risk?
1712
1713  //try and remove
1714  double infty = 1e20;
1715  bool colUpBounded = colUpper < 10000;
1716  bool colLoBounded = colLower > -10000;
1717  bool rowNotLoBounded =  rowLower <= - infty;
1718  bool rowNotUpBounded = rowUpper >= infty;
1719  bool pos =  value > 0;
1720
1721  if(colLoBounded && pos && rowNotUpBounded) {
1722    lb += value * (colsol - colLower);
1723    return 0;
1724  }
1725  else
1726    if(colLoBounded && !pos && rowNotLoBounded) {
1727      ub += value * (colsol - colLower);
1728      return 0;
1729    }
1730    else
1731      if(colUpBounded && !pos && rowNotUpBounded) {
1732        lb += value * (colsol - colUpper);
1733        return 0;
1734      }
1735      else
1736        if(colUpBounded && pos && rowNotLoBounded) {
1737          ub += value * (colsol - colUpper);
1738          return 0;
1739        }
1740  //can not remove coefficient increase it to smallest non zero
1741  if(pos) value = tiny;
1742  else
1743    value = - tiny;
1744  return 1;
1745}
1746
1747/** Get the outer approximation constraints at the point x.
1748*/
1749void
1750OsiTMINLPInterface::getOuterApproximation(OsiCuts &cs, const double * x, 
1751                                          int getObj, const double * x2,
1752                                          double theta, bool global)
1753{
1754  int n,m, nnz_jac_g, nnz_h_lag;
1755  TNLP::IndexStyleEnum index_style;
1756  problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
1757  if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
1758    initializeJacobianArrays();
1759  assert(jRow_ != NULL);
1760  assert(jCol_ != NULL);
1761  vector<double> g(m);
1762  problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
1763  problem_to_optimize_->eval_g(n,x,1,m,g());
1764  //As jacobian is stored by cols fill OsiCuts with cuts
1765  vector<CoinPackedVector>  cuts(nNonLinear_ + 1);
1766  vector<double> lb(nNonLinear_ + 1);
1767  vector<double> ub(nNonLinear_ + 1);
1768
1769  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.
1770  int numCuts = 0;
1771
1772  const double * rowLower = getRowLower();
1773  const double * rowUpper = getRowUpper();
1774  const double * colLower = getColLower();
1775  const double * colUpper = getColUpper();
1776  const double * duals = getRowPrice() + 2 * n;
1777  double infty = getInfinity();
1778  double nlp_infty = infty_;
1779 
1780  for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
1781    if(constTypes_[rowIdx] == TNLP::NON_LINEAR) {
1782      row2cutIdx[rowIdx] = numCuts;
1783      if(rowLower[rowIdx] > - nlp_infty)
1784        lb[numCuts] = rowLower[rowIdx] - g[rowIdx];
1785      else
1786        lb[numCuts] = - infty;
1787      if(rowUpper[rowIdx] < nlp_infty)
1788        ub[numCuts] = rowUpper[rowIdx] - g[rowIdx];
1789      else
1790        ub[numCuts] = infty;
1791      if(rowLower[rowIdx] > -infty && rowUpper[rowIdx] < infty)
1792      {
1793        if(duals[rowIdx] >= 0)// <= inequality
1794          lb[numCuts] = - infty;
1795        if(duals[rowIdx] <= 0)// >= inequality
1796          ub[numCuts] = infty;
1797      }
1798     
1799      numCuts++;
1800    }
1801  }
1802
1803
1804  for(int i = 0 ; i < nnz_jac_g ; i++) {
1805    const int &rowIdx = jRow_[i];
1806    const int & cutIdx = row2cutIdx[ rowIdx ];
1807    if(cutIdx != -1) {
1808      const int &colIdx = jCol_[i];
1809      //"clean" coefficient
1810      if(cleanNnz(jValues_[i],colLower[colIdx], colUpper[colIdx],
1811                  rowLower[rowIdx], rowUpper[rowIdx],
1812                  x[colIdx],
1813                  lb[cutIdx],
1814                  ub[cutIdx], tiny_, veryTiny_)) {
1815        cuts[cutIdx].insert(colIdx,jValues_[i]);
1816        if(lb[cutIdx] > - infty)
1817          lb[cutIdx] += jValues_[i] * x[colIdx];
1818        if(ub[cutIdx] < infty)
1819          ub[cutIdx] += jValues_[i] * x[colIdx];
1820      }
1821    }
1822  }
1823
1824  vector<int> cut2rowIdx(0);
1825  if (IsValid(cutStrengthener_) || oaHandler_->logLevel() > 0) {
1826    cut2rowIdx.resize(numCuts);// Store correspondance between indices of cut and indices of rows. For each cut
1827    for(int rowIdx = 0 ; rowIdx < m ; rowIdx++){
1828       if(row2cutIdx[rowIdx] >= 0){
1829          cut2rowIdx[row2cutIdx[rowIdx]] = rowIdx;
1830       }
1831    }
1832  }
1833
1834  for(int cutIdx = 0; cutIdx < numCuts; cutIdx++) {
1835    //Compute cut violation
1836    if(x2 != NULL) {
1837      double rhs = cuts[cutIdx].dotProduct(x2);
1838      double violation = 0.;
1839      violation = std::max(violation, rhs - ub[cutIdx]);
1840      violation = std::max(violation, lb[cutIdx] - rhs);
1841      if(violation < theta) {
1842        if(oaHandler_->logLevel() > 0)
1843          oaHandler_->message(CUT_NOT_VIOLATED_ENOUGH, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
1844        continue;}
1845        if(oaHandler_->logLevel() > 0)
1846          oaHandler_->message(VIOLATED_OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
1847    }
1848    else if (oaHandler_->logLevel() > 0)
1849      oaHandler_->message(OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<CoinMessageEol;
1850  OsiRowCut newCut;
1851    //    if(lb[i]>-1e20) assert (ub[i]>1e20);
1852
1853    if (IsValid(cutStrengthener_)) {
1854      const int& rowIdx = cut2rowIdx[cutIdx];
1855      bool retval =
1856        cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
1857                                       GetRawPtr(problem_), rowIdx,
1858                                       cuts[cutIdx], lb[cutIdx], ub[cutIdx], g[rowIdx],
1859                                       rowLower[rowIdx], rowUpper[rowIdx],
1860                                       n, x, infty);
1861      if (!retval) {
1862        (*messageHandler()) << "error in cutStrengthener_->ComputeCuts\n";
1863        //exit(-2);
1864      }
1865    }
1866    if(global) {
1867      newCut.setGloballyValidAsInteger(1);
1868    }
1869    //newCut.setEffectiveness(99.99e99);
1870    newCut.setLb(lb[cutIdx]);
1871    newCut.setUb(ub[cutIdx]);
1872    newCut.setRow(cuts[cutIdx]);
1873    if(oaHandler_->logLevel()>2){
1874      oaHandler_->print(newCut);}
1875    cs.insert(newCut);
1876  }
1877
1878  if(getObj == 2 || (getObj && !problem_->hasLinearObjective())) { // Get the objective cuts
1879    vector<double> obj(n);
1880    problem_to_optimize_->eval_grad_f(n, x, 1,obj());
1881    double f;
1882    problem_to_optimize_->eval_f(n, x, 1, f);
1883
1884    CoinPackedVector v;
1885    v.reserve(n);
1886    lb[nNonLinear_] = -f;
1887    ub[nNonLinear_] = -f;
1888    //double minCoeff = 1e50;
1889    for(int i = 0; i<n ; i++) {
1890      if(cleanNnz(obj[i],colLower[i], colUpper[i],
1891          -getInfinity(), 0,
1892          x[i],
1893          lb[nNonLinear_],
1894          ub[nNonLinear_],tiny_, 1e-15)) {
1895        //            minCoeff = min(fabs(obj[i]), minCoeff);
1896        v.insert(i,obj[i]);
1897        lb[nNonLinear_] += obj[i] * x[i];
1898        ub[nNonLinear_] += obj[i] * x[i];
1899      }
1900    }
1901    v.insert(n,-1);
1902    //Compute cut violation
1903    bool genCut = true;
1904    if(x2 != NULL) {
1905      double rhs = v.dotProduct(x2);
1906      double violation = std::max(0., rhs - ub[nNonLinear_]);
1907      if(violation < theta) genCut = false;
1908    }
1909    if(genCut) {
1910      if (IsValid(cutStrengthener_)) {
1911        lb[nNonLinear_] = -infty;
1912        bool retval =
1913          cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
1914                                         GetRawPtr(problem_), -1,
1915                                         v, lb[nNonLinear_], ub[nNonLinear_],
1916                                         ub[nNonLinear_], -infty, 0.,
1917                                         n, x, infty);
1918        if (!retval) {
1919    (*handler_)<< "error in cutStrengthener_->ComputeCuts"<<CoinMessageEol;
1920          //exit(-2);
1921        }
1922      }
1923      OsiRowCut newCut;
1924      if(global)
1925        newCut.setGloballyValidAsInteger(1);
1926      //newCut.setEffectiveness(99.99e99);
1927      newCut.setRow(v);
1928      newCut.setLb(-COIN_DBL_MAX/*Infinity*/);
1929      newCut.setUb(ub[nNonLinear_]);
1930      cs.insert(newCut);
1931    }
1932    }
1933}
1934
1935/** Get a benders cut from solution.*/
1936void
1937OsiTMINLPInterface::getBendersCut(OsiCuts &cs, 
1938                                  bool global){
1939  int n,m, nnz_jac_g, nnz_h_lag;
1940  TNLP::IndexStyleEnum index_style;
1941  problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
1942  if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
1943    initializeJacobianArrays();
1944  assert(jRow_ != NULL);
1945  assert(jCol_ != NULL);
1946  vector<double> g(m);
1947  const double * x = getColSolution();
1948  problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
1949  problem_to_optimize_->eval_g(n,x,1,m,g());
1950  //As jacobian is stored by cols fill OsiCuts with cuts
1951  vector<double> cut(n+1,0.);
1952  vector<bool> keep(m+1,false);
1953  double lb = 0;
1954  double ub = 0;
1955
1956  const double * rowLower = getRowLower();
1957  const double * rowUpper = getRowUpper();
1958  const double * colLower = getColLower();
1959  const double * colUpper = getColUpper();
1960  const double * duals = getRowPrice() + 2 * n;
1961  //double infty = getInfinity();
1962  //double nlp_infty = infty_;
1963 
1964  for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
1965    if(constTypes_[rowIdx] == TNLP::NON_LINEAR && fabs(duals[rowIdx]) > 1e-06)
1966    {
1967      const double & lam = duals[rowIdx];
1968      keep[rowIdx] = true;
1969      assert(lam < 0 || rowUpper[rowIdx] < 1e10);
1970      assert(lam > 0 || rowLower[rowIdx] > -1e10);
1971      if(lam < 0){
1972        assert(rowLower[rowIdx] > -1e10);
1973        ub += lam*(rowLower[rowIdx] -g[rowIdx]);
1974      }
1975      else {
1976        assert(rowUpper[rowIdx] < 1e10);
1977        ub += lam*(rowUpper[rowIdx] -g[rowIdx]);
1978      }
1979    }
1980  }
1981
1982
1983  for(int i = 0 ; i < nnz_jac_g ; i++) {
1984    const int &rowIdx = jRow_[i];
1985    if (!keep[rowIdx]) continue;
1986    const int &colIdx = jCol_[i];
1987    //"clean" coefficient
1988    const double & lam = duals[rowIdx];
1989    double coeff = lam*jValues_[i];
1990    if(cleanNnz(coeff,colLower[colIdx], colUpper[colIdx],
1991          rowLower[rowIdx], rowUpper[rowIdx], x[colIdx], lb,
1992          ub, tiny_, veryTiny_)) {
1993      cut[colIdx] += coeff;
1994      ub += coeff * x[colIdx];
1995    }
1996  }
1997
1998  CoinPackedVector v;
1999  if(!problem_->hasLinearObjective() && 
2000     isProvenOptimal()) { // Get the objective cuts
2001    vector<double> obj(n);
2002    problem_to_optimize_->eval_grad_f(n, x, 1,obj());
2003    double f;
2004    problem_to_optimize_->eval_f(n, x, 1, f);
2005
2006    ub = -f;
2007    //double minCoeff = 1e50;
2008    for(int i = 0; i<n ; i++) {
2009      if(cleanNnz(obj[i],colLower[i], colUpper[i], -getInfinity(), 0,
2010          x[i], lb, ub,tiny_, 1e-15)) {
2011        cut[i] += obj[i];
2012        ub += obj[i] * x[i];
2013      }
2014    }
2015  v.insert(n,-1);
2016  }
2017  for(int i = 0 ; i < n ; i++){
2018    if(fabs(cut[i])>1e-020){
2019      v.insert(i, cut[i]);
2020    }
2021  }
2022  OsiRowCut newCut;
2023  if(global)
2024    newCut.setGloballyValidAsInteger(1);
2025  newCut.setLb(-COIN_DBL_MAX/*Infinity*/);
2026  newCut.setUb(ub);
2027  newCut.setRow(v);
2028  cs.insert(newCut);
2029}
2030
2031/** Get the outer approximation of a single constraint at the point x.
2032*/
2033void
2034OsiTMINLPInterface::getConstraintOuterApproximation(OsiCuts &cs, int rowIdx, 
2035                                                    const double * x, 
2036                                                    const double * x2, bool global)
2037{
2038  double g;
2039  int * indices = new int[getNumCols()];
2040  double * values = new double[getNumCols()];
2041  int nnz;
2042  problem_->eval_grad_gi(getNumCols(), x, 1, rowIdx, nnz, indices, values);
2043  problem_->eval_gi(getNumCols(),x,1, rowIdx, g);
2044
2045  CoinPackedVector cut;
2046  double lb;
2047  double ub;
2048
2049
2050  const double rowLower = getRowLower()[rowIdx];
2051  const double rowUpper = getRowUpper()[rowIdx];
2052  const double * colLower = getColLower();
2053  const double * colUpper = getColUpper();
2054  const double dual = (getRowPrice() + 2 * getNumCols())[rowIdx];
2055  double infty = getInfinity();
2056  double nlp_infty = infty_;
2057 
2058  if(rowLower > - nlp_infty)
2059    lb = rowLower - g;
2060  else
2061    lb = - infty;
2062  if(rowUpper < nlp_infty)
2063    ub = rowUpper - g;
2064  else
2065    ub = infty;
2066  if(rowLower > -infty && rowUpper < infty)
2067  {
2068    if(dual >= 0)// <= inequality
2069      lb = - infty;
2070    if(dual <= 0)// >= inequality
2071      ub = infty;
2072  }
2073
2074  for(int i = 0 ; i < nnz; i++) {
2075     const int &colIdx = indices[i];
2076      //"clean" coefficient
2077      if(cleanNnz(values[i],colLower[colIdx], colUpper[colIdx],
2078                  rowLower, rowUpper,
2079                  x[colIdx],
2080                  lb,
2081                  ub, tiny_, veryTiny_)) {
2082        cut.insert(colIdx,values[i]);
2083        if(lb > - infty)
2084          lb += values[i] * x[colIdx];
2085        if(ub < infty)
2086          ub += values[i] * x[colIdx];
2087    }
2088  }
2089
2090  OsiRowCut newCut;
2091
2092  if(global) {
2093    newCut.setGloballyValidAsInteger(1);
2094  }
2095  //newCut.setEffectiveness(99.99e99);
2096  newCut.setLb(lb);
2097  newCut.setUb(ub);
2098  newCut.setRow(cut);
2099  cs.insert(newCut);
2100
2101  delete [] indices;
2102  delete [] values;
2103}
2104
2105void
2106OsiTMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
2107                                            double a, double s, int L){
2108  if(! IsValid(feasibilityProblem_)) {
2109    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2110  }
2111  feasibilityProblem_->set_use_feasibility_pump_objective(true);
2112  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2113  feasibilityProblem_->setLambda(a);
2114  feasibilityProblem_->setSigma(s);
2115  feasibilityProblem_->setNorm(L);
2116  feasibilityProblem_->set_use_cutoff_constraint(false);
2117  feasibilityProblem_->set_use_local_branching_constraint(false); 
2118  problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
2119  feasibility_mode_ = true;
2120}
2121
2122void
2123OsiTMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
2124                                               double rhs_local_branching_constraint){
2125  if(! IsValid(feasibilityProblem_)) {
2126    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2127  }
2128  feasibilityProblem_->set_use_feasibility_pump_objective(false);
2129  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2130  feasibilityProblem_->set_use_cutoff_constraint(false);
2131  feasibilityProblem_->set_use_local_branching_constraint(true); 
2132  feasibilityProblem_->set_rhs_local_branching_constraint(rhs_local_branching_constraint); 
2133  problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
2134  feasibility_mode_ = true;
2135}
2136
2137void
2138OsiTMINLPInterface::switchToOriginalProblem(){
2139  problem_to_optimize_ = GetRawPtr(problem_);
2140  feasibility_mode_ = false;
2141}
2142
2143double
2144OsiTMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
2145                                            double a, double s, int L)
2146{
2147  if(! IsValid(feasibilityProblem_)) {
2148    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2149  }
2150  feasibilityProblem_->set_use_feasibility_pump_objective(true);
2151  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2152  feasibilityProblem_->setLambda(a);
2153  feasibilityProblem_->setSigma(s);
2154  feasibilityProblem_->setNorm(L);
2155  feasibilityProblem_->set_use_cutoff_constraint(false);
2156  feasibilityProblem_->set_use_local_branching_constraint(false); 
2157  nCallOptimizeTNLP_++;
2158  totalNlpSolveTime_-=CoinCpuTime();
2159  SmartPtr<TNLPSolver> app2 = app_->clone();
2160  app2->options()->SetIntegerValue("print_level", (Index) 0);
2161  optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2162  totalNlpSolveTime_+=CoinCpuTime();
2163  hasBeenOptimized_=true;
2164  return getObjValue();
2165}
2166
2167double
2168OsiTMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
2169                                            int L, double cutoff)
2170{
2171  if(! IsValid(feasibilityProblem_)) {
2172    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2173  }
2174  feasibilityProblem_->set_use_feasibility_pump_objective(true);
2175  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2176  feasibilityProblem_->setLambda(1.0);
2177  feasibilityProblem_->setSigma(0.0);
2178  feasibilityProblem_->setNorm(L);
2179  feasibilityProblem_->set_use_cutoff_constraint(true);
2180  feasibilityProblem_->set_cutoff(cutoff);
2181  feasibilityProblem_->set_use_local_branching_constraint(false); 
2182  nCallOptimizeTNLP_++;
2183  totalNlpSolveTime_-=CoinCpuTime();
2184  SmartPtr<TNLPSolver> app2 = app_->clone();
2185  app2->options()->SetIntegerValue("print_level", (Index) 0);
2186  optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2187  totalNlpSolveTime_+=CoinCpuTime();
2188  hasBeenOptimized_=true;
2189  return getObjValue();
2190}
2191
2192#if 0
2193double
2194OsiTMINLPInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs, bool addOnlyViolated, bool global)
2195{
2196  double ret_val = solveFeasibilityProblem(n, x_bar, inds, 1, 0, 2);
2197  getOuterApproximation(cs, getColSolution(), 0, (addOnlyViolated? x_bar:NULL)
2198                        , global);
2199  return ret_val;
2200}
2201#endif
2202
2203static bool WarnedForNonConvexOa=false;
2204
2205void
2206OsiTMINLPInterface::extractLinearRelaxation(OsiSolverInterface &si, 
2207                                            const double * x, bool getObj
2208                                            )
2209{
2210  int n;
2211  int m;
2212  int nnz_jac_g;
2213  int nnz_h_lag;
2214  TNLP::IndexStyleEnum index_style;
2215  //Get problem information
2216  problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2217
2218  //if not allocated allocate spaced for stroring jacobian
2219  //if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2220    initializeJacobianArrays();
2221
2222  //get Jacobian
2223  problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
2224
2225
2226  vector<double> g(m);
2227  problem_to_optimize_->eval_g(n, x, 1, m, g());
2228
2229  vector<double> rowLow(m);
2230  vector<double> rowUp(m);
2231  vector<int> nonBindings(m);;//store non binding constraints (which are to be removed from OA)
2232  int numNonBindings = 0;
2233  const double * rowLower = getRowLower();
2234  const double * rowUpper = getRowUpper();
2235  const double * colLower = getColLower();
2236  const double * colUpper = getColUpper();
2237  const double * duals = getRowPrice() + 2*n;
2238  assert(m==getNumRows());
2239  double infty = si.getInfinity();
2240  double nlp_infty = infty_;
2241 
2242  for(int i = 0 ; i < m ; i++) {
2243    if(constTypes_[i] == TNLP::NON_LINEAR) {
2244      //If constraint is range not binding prepare to remove it
2245      if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty && fabs(duals[i]) == 0.)
2246      {
2247        nonBindings[numNonBindings++] = i;
2248        continue;
2249      }
2250      else
2251        if(rowLower[i] > - nlp_infty){
2252          rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
2253          if(! WarnedForNonConvexOa && rowUpper[i] < nlp_infty){
2254             messageHandler()->message(WARNING_NON_CONVEX_OA, messages_)<<CoinMessageEol;
2255             WarnedForNonConvexOa = true;
2256          }
2257        }
2258      else
2259        rowLow[i] = - infty;
2260      if(rowUpper[i] < nlp_infty)
2261        rowUp[i] =  (rowUpper[i] - g[i]) + 1e-07;
2262      else
2263        rowUp[i] = infty;
2264     
2265      //If equality or ranged constraint only add one side by looking at sign of dual multiplier
2266      if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty)
2267      {
2268        if(duals[i] >= 0.)// <= inequality
2269          rowLow[i] = - infty;
2270        if(duals[i] <= 0.)// >= inequality
2271          rowUp[i] = infty;
2272      }
2273    }
2274    else {
2275      if(rowLower[i] > -nlp_infty){
2276      //   printf("Lower %g ", rowLower[i]);
2277         rowLow[i] = (rowLower[i] - g[i]);
2278      }
2279      else
2280        rowLow[i] = - infty;
2281      if(rowUpper[i] < nlp_infty){
2282      //   printf("Upper %g ", rowUpper[i]);
2283         rowUp[i] =  (rowUpper[i] - g[i]);
2284      }
2285      else
2286        rowUp[i] = infty;
2287    }
2288  }
2289
2290 
2291 
2292  //Then convert everything to a CoinPackedMatrix
2293  //Go through values, clean coefficients and fix bounds
2294  for(int i = 0 ; i < nnz_jac_g ; i++) {
2295    if(constTypes_[jRow_[i]] != TNLP::LINEAR){//For linear just copy is fine.
2296       if(//For other clean tinys
2297       cleanNnz(jValues_[i],colLower[jCol_[i]], colUpper[jCol_[i]],
2298                rowLower[jRow_[i]], rowUpper[jRow_[i]],
2299                x[jCol_[i]],
2300                rowLow[jRow_[i]],
2301                rowUp[jRow_[i]], tiny_, veryTiny_)) {     
2302          rowLow[jRow_[i]] += jValues_[i] * x[jCol_ [i]];
2303          rowUp[jRow_[i]] += jValues_[i] *x[jCol_[i]];
2304       }
2305    }
2306    else {
2307      double value = jValues_[i] * getColSolution()[jCol_[i]];
2308      rowLow[jRow_[i]] += value;
2309      rowUp[jRow_[i]] += value;
2310    } 
2311  }
2312  CoinPackedMatrix mat(true, jRow_, jCol_, jValues_, nnz_jac_g);
2313  mat.setDimensions(m,n); // In case matrix was empty, this should be enough
2314 
2315  //remove non-bindings equality constraints
2316  mat.deleteRows(numNonBindings, nonBindings());
2317
2318  int numcols=getNumCols();
2319  vector<double> obj(numcols);
2320  for(int i = 0 ; i < numcols ; i++)
2321    obj[i] = 0.;
2322 
2323 
2324  si.loadProblem(mat, getColLower(), getColUpper(), obj(), rowLow(), rowUp());
2325  for(int i = 0 ; i < getNumCols() ; i++) {
2326    if(isInteger(i))
2327      si.setInteger(i);
2328  }
2329  if(getObj) {
2330     bool addObjVar = false;
2331     if(problem_->hasLinearObjective()){
2332       double zero;
2333       vector<double> x0(n,0.);
2334       problem_to_optimize_->eval_f(n, x0(), 1, zero);
2335       si.setDblParam(OsiObjOffset, -zero);
2336       //Copy the linear objective and don't create a dummy variable.
2337       problem_to_optimize_->eval_grad_f(n, x, 1,obj());
2338       si.setObjective(obj());
2339    }
2340    else {
2341      addObjVar = true;
2342   }
2343   
2344   if(addObjVar){
2345      addObjectiveFunction(si, x);
2346    }
2347  }
2348//  si.writeMpsNative("OA.mps",NULL, NULL, 1);
2349}
2350
2351void 
2352OsiTMINLPInterface::addObjectiveFunction(OsiSolverInterface &si, 
2353                                         const double *x){
2354  const double * colLower = getColLower();
2355  const double * colUpper = getColUpper();
2356  int numcols = getNumCols();
2357  assert(numcols == si.getNumCols() );
2358  vector<double> obj(numcols);
2359      problem_to_optimize_->eval_grad_f(numcols, x, 1,obj());
2360      //add variable alpha
2361      //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
2362      CoinPackedVector a;
2363      si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
2364 
2365      // Now get the objective cuts
2366      // get the gradient, pack it and add the cut
2367      double ub;
2368      problem_to_optimize_->eval_f(numcols, x, 1, ub);
2369      ub*=-1;
2370      double lb = -1e300;
2371      CoinPackedVector objCut;
2372      CoinPackedVector * v = &objCut;
2373      v->reserve(numcols+1);
2374      for(int i = 0; i<numcols ; i++) {
2375       if(si.getNumRows())
2376       {
2377        if(cleanNnz(obj[i],colLower[i], colUpper[i],
2378            -getInfinity(), 0,
2379            x[i],
2380            lb,
2381            ub, tiny_, veryTiny_)) {
2382          v->insert(i,obj[i]);
2383          lb += obj[i] * x[i];
2384          ub += obj[i] * x[i];
2385        }
2386       }
2387       else //Unconstrained problem can not put clean coefficient
2388       {
2389           if(cleanNnz(obj[i],colLower[i], colUpper[i],
2390            -getInfinity(), 0,
2391            x[i],
2392            lb,
2393            ub, 1e-03, 1e-08)) {
2394          v->insert(i,obj[i]);
2395          lb += obj[i] * x[i];
2396          ub += obj[i] * x[i];
2397           }
2398       }
2399      }
2400    v->insert(numcols,-1);
2401    si.addRow(objCut, lb, ub);
2402    }
2403
2404/** Add a collection of linear cuts to problem formulation.*/
2405void 
2406OsiTMINLPInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
2407{ 
2408  if(numberCuts)
2409    freeCachedRowRim();
2410  const OsiRowCut ** cutsPtrs = new const OsiRowCut*[numberCuts];
2411  for(int i = 0 ; i < numberCuts ; i++)
2412  {
2413    cutsPtrs[i] = &cuts[i];
2414  }
2415  problem_->addCuts(numberCuts, cutsPtrs);
2416  delete [] cutsPtrs;
2417}
2418
2419void
2420OsiTMINLPInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
2421    const char * whereFrom)
2422{
2423  totalNlpSolveTime_-=CoinCpuTime();
2424  if(warmStarted)
2425    optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
2426  else
2427    optimizationStatus_ = app_->OptimizeTNLP(GetRawPtr(problem_to_optimize_));
2428  totalNlpSolveTime_+=CoinCpuTime();
2429  nCallOptimizeTNLP_++;
2430  hasBeenOptimized_ = true;
2431 
2432 
2433  //Options should have been printed if not done already turn off Ipopt output
2434  if(!hasPrintedOptions) {
2435    hasPrintedOptions = 1;
2436    //app_->Options()->SetIntegerValue("print_level",0, true, true);
2437    app_->options()->SetStringValue("print_user_options","no", false, true);
2438  }
2439 
2440  bool otherDisagree = false ;
2441#if 0
2442  if(optimizationStatus_ == TNLPSolver::notEnoughFreedom)//Too few degrees of freedom
2443  {
2444    (*messageHandler())<<"Too few degrees of freedom...."<<CoinMessageEol;
2445    int numrows = getNumRows();
2446    int numcols = getNumCols();
2447   
2448    const double * colLower = getColLower();
2449    const double * colUpper = getColUpper();
2450   
2451   
2452    const double * rowLower = getRowLower();
2453    const double * rowUpper = getRowUpper();
2454   
2455    int numberFixed = 0;
2456    for(int i = 0 ; i < numcols ; i++)
2457    {
2458      if(colUpper[i] - colLower[i] <= INT_BIAS)
2459            {
2460              numberFixed++;
2461            }
2462    }
2463    int numberEqualities = 0;
2464    for(int i = 0 ; i < numrows ; i++)
2465    {
2466      if(rowUpper[i] - rowLower[i] <= INT_BIAS)
2467            {
2468              numberEqualities++;
2469            }     
2470    }
2471    if(numcols - numberFixed > numberEqualities || numcols < numberEqualities)
2472    {
2473      std::string probName;
2474      getStrParam(OsiProbName, probName);
2475      throw newUnsolvedError(app_->errorCode(), problem_, probName);
2476    }
2477    double * saveColLow = CoinCopyOfArray(getColLower(), getNumCols());
2478    double * saveColUp = CoinCopyOfArray(getColUpper(), getNumCols());
2479   
2480    for(int i = 0; i < numcols && numcols - numberFixed <= numberEqualities ; i++)
2481    {
2482      if(colUpper[i] - colLower[i] <= INT_BIAS)
2483            {
2484              setColLower(i, saveColLow[i]-1e-06);
2485              setColUpper(i, saveColLow[i]+1e-06);
2486              numberFixed--;
2487            }
2488    }
2489    solveAndCheckErrors(warmStarted, throwOnFailure, whereFrom);
2490    //restore
2491    for(int i = 0; i < numcols && numcols - numberFixed < numrows ; i++)
2492    {
2493      problem_->SetVariableLowerBound(i,saveColLow[i]);
2494      problem_->SetVariableUpperBound(i,saveColUp[i]);
2495    }
2496    delete [] saveColLow;
2497    delete [] saveColUp;
2498    return;
2499  }
2500  else
2501#endif
2502    if(!app_->isRecoverable(optimizationStatus_))//Solver failed and the error can not be recovered, throw it
2503    {
2504      std::string probName;
2505      getStrParam(OsiProbName, probName);
2506      throw newUnsolvedError(app_->errorCode(), problem_, probName);
2507    }
2508    else if(testOthers_ && !app_->isError(optimizationStatus_)){
2509      Ipopt::SmartPtr<TMINLP2TNLP> problem_copy = problem_->clone();
2510      //Try other solvers and see if they agree
2511      int f =1;
2512      for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2513          i != debug_apps_.end() ; i++){
2514        TNLPSolver::ReturnStatus otherStatus = (*i)->OptimizeTNLP(GetRawPtr(problem_copy));
2515       messageHandler()->message(LOG_LINE, messages_)
2516         <<'d'<<f++<<statusAsString(otherStatus)<<problem_copy->obj_value()
2517         <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
2518        if(!(*i)->isError(otherStatus)){
2519           CoinRelFltEq eq(1e-05);
2520           if(otherStatus != optimizationStatus_){
2521             otherDisagree = true;
2522             messageHandler()->message(SOLVER_DISAGREE_STATUS, messages_)
2523             <<app_->solverName()<<statusAsString()
2524             <<(*i)->solverName()<<statusAsString(otherStatus)<<CoinMessageEol; 
2525           }
2526           else if(isProvenOptimal() && !eq(problem_->obj_value(),problem_copy->obj_value()))
2527           {
2528             otherDisagree = true;
2529             messageHandler()->message(SOLVER_DISAGREE_VALUE, messages_)
2530             <<app_->solverName()<<problem_->obj_value()
2531             <<(*i)->solverName()<<problem_copy->obj_value()<<CoinMessageEol; 
2532           }
2533        }
2534     }
2535  }
2536  else if(app_->isError(optimizationStatus_) && ! debug_apps_.empty()){
2537      int f =1;
2538      for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2539          i != debug_apps_.end() && app_->isError(optimizationStatus_) ; i++){
2540        optimizationStatus_ = (*i)->OptimizeTNLP(GetRawPtr(problem_));
2541        messageHandler()->message(LOG_LINE, messages_)
2542          <<'d'<<f++<<statusAsString(optimizationStatus_)<<problem_->obj_value()
2543          <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
2544      }
2545  }
2546  try{
2547    totalIterations_ += app_->IterationCount();
2548  }
2549  catch(SimpleError &E)
2550  {
2551    if (throwOnFailure)//something failed throw
2552    {
2553      throw SimpleError("No statistics available from Ipopt",whereFrom);
2554    }
2555    else {
2556      return;
2557    }
2558  }
2559  if(problem_->hasUpperBoundingObjective()){//Check if solution is integer and recompute objective value using alternative objective function
2560    const double * sol = getColSolution();
2561    bool integerSol = true;
2562    double intTol = 1e-08;
2563    if(objects()){
2564      int nObjects = numberObjects();
2565      OsiObject ** object = objects();
2566      for(int i = 0 ; i< nObjects ; i++){
2567        int dummy;
2568        if(object[i]->infeasibility(this,dummy) > intTol)
2569        {
2570          integerSol=false;
2571          break;
2572        }
2573      }
2574    }
2575    else{//Only works for integer constraints
2576      int numcols = getNumCols();
2577      for(int i = 0 ; i < numcols ; i++){
2578        if(isInteger(i) || isBinary(i)){
2579          if(fabs(sol[i] - floor(sol[i]+0.5)) > intTol){
2580            integerSol = false;
2581            break;
2582          }
2583        }
2584      }
2585    }
2586    if(integerSol&&isProvenOptimal()){
2587      double help= problem_->evaluateUpperBoundingFunction(sol);
2588     
2589
2590      OsiAuxInfo * auxInfo = getAuxiliaryInfo();
2591      Bonmin::AuxInfo * bonInfo = dynamic_cast<Bonmin::AuxInfo *>(auxInfo);
2592      if(bonInfo!=0)
2593      {
2594       
2595        if(help<bonInfo->bestObj2())
2596        {
2597          bonInfo->setBestObj2(help);
2598          bonInfo->setBestSolution2(getNumCols(), const_cast<double *>(getColSolution()));
2599
2600           messageHandler()->message(ALTERNATE_OBJECTIVE, messages_)
2601           <<help<<CoinMessageEol;
2602        }
2603      }
2604      else {
2605        printf("\nWARNING: the algorithm selected does not consider the second objective function\n");
2606      }
2607    }
2608  }
2609  messageHandler()->message(IPOPT_SUMMARY, messages_)
2610    <<whereFrom<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
2611 
2612  if((nCallOptimizeTNLP_ % 20) == 1)
2613    messageHandler()->message(LOG_HEAD, messages_)<<CoinMessageEol;
2614 
2615 
2616  if ( (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) ||
2617       ( otherDisagree )){
2618    messageHandler()->message(SUSPECT_PROBLEM,
2619                              messages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
2620    std::string subProbName;
2621    getStrParam(OsiProbName, subProbName);
2622    std::ostringstream os;
2623    os<<"_"<<nCallOptimizeTNLP_;
2624    subProbName+=os.str();
2625    problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
2626  }
2627 
2628}
2629
2630////////////////////////////////////////////////////////////////////
2631// Solve Methods                                                  //
2632////////////////////////////////////////////////////////////////////
2633void OsiTMINLPInterface::initialSolve()
2634{
2635   initialSolve("");
2636}
2637
2638////////////////////////////////////////////////////////////////////
2639// Solve Methods                                                  //
2640////////////////////////////////////////////////////////////////////
2641void OsiTMINLPInterface::initialSolve(const char * whereFrom)
2642{
2643  assert(IsValid(app_));
2644  assert(IsValid(problem_));
2645
2646  // Discard warmstart_ if we had one
2647  delete warmstart_;
2648  warmstart_ = NULL;
2649 
2650  if(!hasPrintedOptions) {
2651    int printOptions;
2652    app_->options()->GetEnumValue("print_user_options",printOptions,app_->prefix());
2653    if(printOptions)
2654      app_->options()->SetStringValue("print_user_options","yes",true,true);
2655  }
2656  if(exposeWarmStart_)
2657    app_->disableWarmStart(); 
2658  solveAndCheckErrors(0,1,"initialSolve");
2659 
2660  //Options should have been printed if not done already turn off Ipopt output
2661  if(!hasPrintedOptions) {
2662    hasPrintedOptions = 1;
2663    app_->options()->SetStringValue("print_user_options","no");
2664    app_->options()->SetIntegerValue("print_level",0);
2665  }
2666 
2667  messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
2668                                                      <<statusAsString()
2669                                                      <<getObjValue()
2670                                                      <<app_->IterationCount()
2671                                                      <<app_->CPUTime()
2672                                                      <<whereFrom<<CoinMessageEol;
2673 
2674  int numRetry = firstSolve_ ? numRetryInitial_ : numRetryResolve_;
2675  if(isAbandoned()) {
2676    resolveForRobustness(numRetryUnsolved_);
2677  }
2678  else if(numRetry)
2679    {
2680      resolveForCost(numRetry, numRetryInitial_ > 0);
2681      /** Only do it once at the root.*/
2682      numRetryInitial_ = 0;
2683    }
2684  firstSolve_ = false;
2685
2686  // if warmstart_ is not empty then had to use resolveFor... and that created
2687  // the warmstart at the end, and we have nothing to do here. Otherwise...
2688  if (! warmstart_ && ! isAbandoned()) {
2689    if (exposeWarmStart_) {
2690      warmstart_ = app_->getWarmStart(problem_);
2691    }
2692  }
2693}
2694
2695/** Resolve the continuous relaxation after problem modification. */
2696void
2697OsiTMINLPInterface::resolve(){
2698   resolve("");
2699}
2700/** Resolve the continuous relaxation after problem modification.
2701 * \note for Ipopt, same as resolve */
2702void
2703OsiTMINLPInterface::resolve(const char * whereFrom)
2704{
2705  assert(IsValid(app_));
2706  assert(IsValid(problem_));
2707  int has_warmstart = warmstart_ == NULL ? 0 : 1;
2708  if(warmstart_ == NULL) has_warmstart = 0;
2709  else if(!app_->warmStartIsValid(warmstart_)) has_warmstart = 1;
2710  else has_warmstart = 2;
2711  if (has_warmstart < 2) {
2712    // empty (or unrecognized) warmstart
2713    initialSolve(whereFrom);
2714    return;
2715  }
2716  app_->setWarmStart(warmstart_, problem_);
2717  delete warmstart_;
2718  warmstart_ = NULL;
2719
2720  if (INT_BIAS > 0.) {
2721    app_->options()->SetStringValue("warm_start_same_structure", "yes");
2722  }
2723  else {
2724    app_->options()->SetStringValue("warm_start_same_structure", "no");
2725  }
2726
2727  if(problem_->duals_init() != NULL)
2728    app_->enableWarmStart();
2729  else app_->disableWarmStart();
2730  solveAndCheckErrors(1,1,"resolve");
2731 
2732  messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
2733                                                <<statusAsString()
2734                                                <<getObjValue()
2735                                                <<app_->IterationCount()
2736                                                <<app_->CPUTime()
2737                                                <<whereFrom
2738                                                <<"totot"
2739                                                <<CoinMessageEol;
2740 
2741  if(isAbandoned()) {
2742    resolveForRobustness(numRetryUnsolved_);
2743  }
2744  else if(numRetryResolve_ ||
2745          (numRetryInfeasibles_ && isProvenPrimalInfeasible() ))
2746    resolveForCost(std::max(numRetryResolve_, numRetryInfeasibles_), 0);
2747
2748  // if warmstart_ is not empty then had to use resolveFor... and that created
2749  // the warmstart at the end, and we have nothing to do here. Otherwise...
2750  if (! warmstart_ && ! isAbandoned()) {
2751    if (exposeWarmStart_) {
2752      warmstart_ = app_->getWarmStart(problem_);
2753    }
2754  }
2755}
2756
2757
2758////////////////////////////////////////////////////////////////
2759// Methods returning info on how the solution process terminated  //
2760////////////////////////////////////////////////////////////////
2761/// Are there a numerical difficulties?
2762bool OsiTMINLPInterface::isAbandoned() const
2763{
2764  return (
2765        (optimizationStatus_==TNLPSolver::iterationLimit)||
2766        (optimizationStatus_==TNLPSolver::computationError)||
2767        (optimizationStatus_==TNLPSolver::illDefinedProblem)||
2768        (optimizationStatus_==TNLPSolver::illegalOption)||
2769        (optimizationStatus_==TNLPSolver::externalException)|
2770        (optimizationStatus_==TNLPSolver::exception)
2771      );
2772}
2773
2774/// Is optimality proven?
2775bool OsiTMINLPInterface::isProvenOptimal() const
2776{
2777  return (optimizationStatus_==TNLPSolver::solvedOptimal) ||
2778          (optimizationStatus_==TNLPSolver::solvedOptimalTol);
2779}
2780/// Is primal infeasiblity proven?
2781bool OsiTMINLPInterface::isProvenPrimalInfeasible() const
2782{
2783  return (optimizationStatus_ == TNLPSolver::provenInfeasible);
2784}
2785/// Is dual infeasiblity proven?
2786bool OsiTMINLPInterface::isProvenDualInfeasible() const
2787{
2788  return (optimizationStatus_ == TNLPSolver::unbounded);
2789}
2790/// Is the given primal objective limit reached?
2791bool OsiTMINLPInterface::isPrimalObjectiveLimitReached() const
2792{
2793  (*handler_)<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2794  return 0;
2795}
2796/// Is the given dual objective limit reached?
2797bool OsiTMINLPInterface::isDualObjectiveLimitReached() const
2798{
2799  //  (*messageHandler_)<<"Warning : isDualObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2800  return (optimizationStatus_==TNLPSolver::unbounded);
2801
2802}
2803/// Iteration limit reached?
2804bool OsiTMINLPInterface::isIterationLimitReached() const
2805{
2806  return (optimizationStatus_==TNLPSolver::iterationLimit);
2807}
2808
2809void
2810OsiTMINLPInterface::extractInterfaceParams()
2811{
2812  if (IsValid(app_)) {
2813    int logLevel;
2814    app_->options()->GetIntegerValue("nlp_log_level", logLevel,app_->prefix());
2815    messageHandler()->setLogLevel(logLevel);
2816
2817#ifdef COIN_HAS_FILTERSQP
2818    FilterSolver * filter = dynamic_cast<FilterSolver *>(GetRawPtr(app_));
2819
2820    bool is_given =
2821#endif
2822      app_->options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,app_->prefix());
2823
2824#ifdef COIN_HAS_FILTERSQP
2825    if(filter && !is_given){
2826      // overwriting default for filter
2827      maxRandomRadius_ = 10.;
2828    }
2829#endif
2830   
2831   int oaCgLogLevel = 0;
2832   app_->options()->GetIntegerValue("oa_cuts_log_level", oaCgLogLevel,app_->prefix());
2833   oaHandler_->setLogLevel(oaCgLogLevel); 
2834   
2835    int exposeWs = false;
2836    app_->options()->GetEnumValue("warm_start", exposeWs, app_->prefix());
2837    setExposeWarmStart(exposeWs > 0);
2838   
2839    app_->options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,app_->prefix());
2840    app_->options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,app_->prefix());
2841    app_->options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,app_->prefix());
2842    app_->options()->GetIntegerValue("num_resolve_at_infeasibles", numRetryInfeasibles_,app_->prefix());
2843    app_->options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,app_->prefix());
2844    app_->options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,app_->prefix());
2845    app_->options()->GetNumericValue
2846    ("warm_start_bound_frac" ,pushValue_,app_->prefix());
2847    app_->options()->GetNumericValue("tiny_element",tiny_,app_->prefix());
2848    app_->options()->GetNumericValue("very_tiny_element",veryTiny_,app_->prefix());
2849    app_->options()->GetNumericValue("random_point_perturbation_interval",max_perturbation_,app_->prefix());
2850    app_->options()->GetEnumValue("random_point_type",randomGenerationType_,app_->prefix());
2851    int cut_strengthening_type;
2852    app_->options()->GetEnumValue("cut_strengthening_type", cut_strengthening_type,app_->prefix());
2853
2854    if (cut_strengthening_type != CS_None) {
2855      // TNLP solver to be used in the cut strengthener
2856      cutStrengthener_ = new CutStrengthener(app_->clone(), app_->options());
2857    }
2858  }
2859}
2860
2861void
2862OsiTMINLPInterface::SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)
2863{
2864  strong_branching_solver_ = strong_branching_solver;
2865}
2866
2867  //#define STRONG_COMPARE
2868#ifdef STRONG_COMPARE
2869  static double objorig;
2870#endif
2871
2872void
2873OsiTMINLPInterface::markHotStart()
2874{
2875  if (IsValid(strong_branching_solver_)) {
2876#ifdef STRONG_COMPARE
2877    // AWDEBUG
2878    OsiSolverInterface::markHotStart();
2879    objorig = getObjValue();
2880#endif
2881    optimizationStatusBeforeHotStart_ = optimizationStatus_;
2882    strong_branching_solver_->markHotStart(this);
2883  }
2884  else {
2885    // Default Implementation
2886    OsiSolverInterface::markHotStart();
2887  }
2888}
2889
2890void
2891OsiTMINLPInterface::solveFromHotStart()
2892{
2893  if (IsValid(strong_branching_solver_)) {
2894#ifdef STRONG_COMPARE
2895    // AWDEBUG
2896    OsiSolverInterface::solveFromHotStart();
2897    double obj_nlp = getObjValue() - objorig;
2898#endif
2899    optimizationStatus_ = strong_branching_solver_->solveFromHotStart(this);
2900    hasBeenOptimized_ = true;
2901#ifdef STRONG_COMPARE
2902    double obj_other = getObjValue() - objorig;
2903    printf("AWDEBUG: Strong Branching results: NLP = %15.8e Other = %15.8e\n",
2904           obj_nlp, obj_other);
2905#endif
2906  }
2907  else {
2908    // Default Implementation
2909    OsiSolverInterface::solveFromHotStart();
2910  }
2911}
2912
2913void
2914OsiTMINLPInterface::unmarkHotStart()
2915{
2916  if (IsValid(strong_branching_solver_)) {
2917#ifdef STRONG_COMPARE
2918    // AWDEBUG
2919    OsiSolverInterface::unmarkHotStart();
2920#endif
2921    strong_branching_solver_->unmarkHotStart(this);
2922    optimizationStatus_ = optimizationStatusBeforeHotStart_;
2923  }
2924  else {
2925    // Default Implementation
2926    OsiSolverInterface::unmarkHotStart();
2927  }
2928}
2929
2930const double * OsiTMINLPInterface::getObjCoefficients() const
2931{
2932  const int n = getNumCols();
2933  delete [] obj_;
2934  obj_ = NULL;
2935  obj_ = new double[n];
2936
2937  bool new_x = true;
2938  const double* x_sol = problem_->x_sol();
2939  bool retval = problem_->eval_grad_f(n, x_sol, new_x, obj_);
2940 
2941  if (!retval) {
2942    // Let's see if that happens - it will cause a crash
2943    fprintf(stderr, "ERROR WHILE EVALUATING GRAD_F in OsiTMINLPInterface::getObjCoefficients()\n");
2944    delete [] obj_;
2945    obj_ = NULL;
2946  }
2947
2948  return obj_;
2949}
2950
2951  /** Sets the TMINLP2TNLP to be used by the interface.*/
2952  void 
2953  OsiTMINLPInterface::use(Ipopt::SmartPtr<TMINLP2TNLP> tminlp2tnlp){
2954     problem_ = tminlp2tnlp;
2955     feasibilityProblem_->use(GetRawPtr(tminlp2tnlp));}
2956
2957}/** end namespace Bonmin*/
2958
Note: See TracBrowser for help on using the repository browser.