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

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

More fixes for unconstrained problems with windows

File size: 87.7 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                             "");
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");
273  ADD_MSG(LOG_FIRST_LINE, std_m, 1, 
274          "    %-8d %-11s %-23.16g %-8d %-3g");
275  ADD_MSG(LOG_LINE, std_m, 1, " %c  r%-7d %-11s %-23.16g %-8d %-3g");
276
277  ADD_MSG(ALTERNATE_OBJECTIVE, std_m, 1, "Objective value recomputed with alternate objective: %g.");
278 
279  ADD_MSG(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, warn_m, 1, 
280         "resolve called before any call to initialSol  can not use warm starts.");
281  ADD_MSG(ERROR_NO_TNLPSOLVER, warn_m, 1,"Can not parse options when no IpApplication has been created");
282  ADD_MSG(WARNING_NON_CONVEX_OA, warn_m, 1,
283          "OA on non-convex constraint is very experimental.");                         
284  ADD_MSG(SOLVER_DISAGREE_STATUS, warn_m, 1, "%s says problem %s, %s says %s.");
285  ADD_MSG(SOLVER_DISAGREE_VALUE, warn_m, 1, "%s gives objective %.16g, %s gives %.16g.");
286
287}
288
289
290void 
291OsiTMINLPInterface::OaMessageHandler::print(OsiRowCut &row){
292  FILE * fp = filePointer();
293  const int &n = row.row().getNumElements();
294  fprintf(fp,"Row cut has %d elements. Lower bound: %g, upper bound %g.\n", n, row.lb(), row.ub());
295  const int * idx = row.row().getIndices();
296  const double * val = row.row().getElements();
297  for(int i = 0 ; i < n ; i++){
298    fprintf(fp,"%g, x%d",val[i], idx[i]);
299    if(i && i % 7 == 0)
300      fprintf(fp,"\n");
301  } 
302}
303
304OsiTMINLPInterface::OaMessages::OaMessages(): CoinMessages((int) OA_MESSAGES_DUMMY_END){
305   strcpy(source_,"OaCg");
306   ADD_MSG(VIOLATED_OA_CUT_GENERATED, std_m, 1,"Row %d, cut violation is %g: Outer approximation cut generated.");
307
308   ADD_MSG(CUT_NOT_VIOLATED_ENOUGH, std_m, 2,"Row %d, cut violation is %g: Outer approximation cut not generated.");
309
310   ADD_MSG(OA_CUT_GENERATED, std_m, 1,"Row %d: Outer approximation cut not generated.");
311}
312bool OsiTMINLPInterface::hasPrintedOptions=0;
313
314////////////////////////////////////////////////////////////////////
315// Constructors and desctructors                                  //
316////////////////////////////////////////////////////////////////////
317/// Default Constructor
318OsiTMINLPInterface::OsiTMINLPInterface():
319    OsiSolverInterface(),
320    tminlp_(NULL),
321    problem_(NULL),
322    problem_to_optimize_(NULL),
323    feasibility_mode_(false),
324    app_(NULL),
325    debug_apps_(),
326    testOthers_(false),
327    warmstart_(NULL),
328    rowsense_(NULL),
329    rhs_(NULL),
330    rowrange_(NULL),
331    reducedCosts_(NULL),
332    OsiDualObjectiveLimit_(1e200),
333    hasVarNamesFile_(true),
334    nCallOptimizeTNLP_(0),
335    totalNlpSolveTime_(0),
336    totalIterations_(0),
337    maxRandomRadius_(1e08),
338    randomGenerationType_(0),
339    max_perturbation_(COIN_DBL_MAX),
340    pushValue_(1e-02),
341    numRetryInitial_(-1),
342    numRetryResolve_(-1),
343    numRetryInfeasibles_(-1),
344    numRetryUnsolved_(1),
345    messages_(),
346    pretendFailIsInfeasible_(0),
347    hasContinuedAfterNlpFailure_(false),
348    numIterationSuspect_(-1),
349    hasBeenOptimized_(false),
350    obj_(NULL),
351    feasibilityProblem_(NULL),
352    jRow_(NULL),
353    jCol_(NULL),
354    jValues_(NULL),
355    nnz_jac(0),
356    constTypes_(NULL),
357    nNonLinear_(0),
358    tiny_(1e-8),
359    veryTiny_(1e-20),
360    infty_(1e100),
361    exposeWarmStart_(false),
362    firstSolve_(true),
363    cutStrengthener_(NULL),
364    oaMessages_(),
365    oaHandler_(NULL)
366{
367   oaHandler_ = new OaMessageHandler;
368   oaHandler_->setLogLevel(0);
369}
370
371void 
372OsiTMINLPInterface::initialize(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
373                               Ipopt::SmartPtr<Ipopt::OptionsList> options,
374                               Ipopt::SmartPtr<Ipopt::Journalist> journalist,
375                               const std::string & prefix,
376                               Ipopt::SmartPtr<TMINLP> tminlp){
377  if(!IsValid(app_))
378     createApplication(roptions, options, journalist, prefix);
379  setModel(tminlp); 
380}
381
382void OsiTMINLPInterface::setSolver(Ipopt::SmartPtr<TNLPSolver> app){
383  app_ = app;
384  }
385
386
387void
388OsiTMINLPInterface::createApplication(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
389                                      Ipopt::SmartPtr<Ipopt::OptionsList> options,
390                                      Ipopt::SmartPtr<Ipopt::Journalist> journalist,
391                                      const std::string & prefix
392                                      )
393{
394  assert(!IsValid(app_));
395  int ival;
396  options->GetEnumValue("nlp_solver", ival, prefix.c_str());
397  Solver s = (Solver) ival;
398  if(s == EFilterSQP){
399    testOthers_ = false;;
400#ifdef COIN_HAS_FILTERSQP
401    app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
402#else
403   throw SimpleError("createApplication",
404                     "Bonmin not configured to run with FilterSQP.");
405#endif   
406
407#ifdef COIN_HAS_IPOPT
408   debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix)); 
409#endif
410  }
411  else if(s == EIpopt){
412    testOthers_ = false;
413#ifdef COIN_HAS_IPOPT
414    app_ = new IpoptSolver(roptions, options, journalist, prefix);
415#else
416   throw SimpleError("createApplication",
417                     "Bonmin not configured to run with Ipopt.");
418#endif
419#ifdef COIN_HAS_FILTERSQP
420    debug_apps_.push_back(new Bonmin::FilterSolver(roptions, options, journalist, prefix));
421#endif
422  }
423  else if(s == EAll){
424#ifdef COIN_HAS_FILTERSQP
425    app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
426#else
427   throw SimpleError("createApplication",
428                     "Bonmin not configured to run with Ipopt.");
429#endif
430#ifdef COIN_HAS_IPOPT
431   debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix)); 
432#endif
433    testOthers_ = true;
434  }
435  if (!app_->Initialize("")) {
436    throw CoinError("Error during initialization of app_","createApplication", "OsiTMINLPInterface");
437  }
438  for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin() ; 
439      i != debug_apps_.end() ; i++){
440    (*i)->Initialize("");
441  }
442  extractInterfaceParams();
443 
444}
445
446/** Facilitator to allocate a tminlp and an application. */
447void
448OsiTMINLPInterface::setModel(SmartPtr<TMINLP> tminlp)
449{
450  assert(IsValid(tminlp));
451  tminlp_ = tminlp;
452  problem_ = new TMINLP2TNLP(tminlp_);
453  feasibilityProblem_ = new TNLP2FPNLP
454        (SmartPtr<TNLP>(GetRawPtr(problem_)));
455  if(feasibility_mode_){
456    problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
457  }
458  else {
459    problem_to_optimize_ = GetRawPtr(problem_);
460  }
461}
462
463
464
465void
466OsiTMINLPInterface::readOptionFile(const std::string & fileName)
467{
468  if(IsValid(app_)){
469  std::ifstream is;
470  if (fileName != "") {
471    try {
472      is.open(fileName.c_str());
473    }
474    catch(std::bad_alloc) {
475      throw CoinError("Not enough memory to open option file.\n", "readOptionFile", "OsiTMINLPInterface");
476    }
477  }
478  options()->ReadFromStream(*app_->journalist(), is);
479  extractInterfaceParams();
480  }
481}
482
483/// Copy constructor
484OsiTMINLPInterface::OsiTMINLPInterface (const OsiTMINLPInterface &source):
485    OsiSolverInterface(source),
486    tminlp_(source.tminlp_),
487    problem_(NULL),
488    problem_to_optimize_(NULL),
489    feasibility_mode_(source.feasibility_mode_),
490    rowsense_(NULL),
491    rhs_(NULL),
492    rowrange_(NULL),
493    reducedCosts_(NULL),
494    OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
495    hasVarNamesFile_(source.hasVarNamesFile_),
496    nCallOptimizeTNLP_(0),
497    totalNlpSolveTime_(0),
498    totalIterations_(0),
499    maxRandomRadius_(source.maxRandomRadius_),
500    randomGenerationType_(source.randomGenerationType_),
501    max_perturbation_(source.max_perturbation_),
502    pushValue_(source.pushValue_),
503    numRetryInitial_(source.numRetryInitial_),
504    numRetryResolve_(source.numRetryResolve_),
505    numRetryInfeasibles_(source.numRetryInfeasibles_),
506    numRetryUnsolved_(source.numRetryUnsolved_),
507    messages_(),
508    pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
509    hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
510    numIterationSuspect_(source.numIterationSuspect_),
511    hasBeenOptimized_(source.hasBeenOptimized_),
512    obj_(NULL),
513    feasibilityProblem_(NULL),
514    jRow_(NULL),
515    jCol_(NULL),
516    jValues_(NULL),
517    nnz_jac(source.nnz_jac),
518    constTypes_(NULL),
519//    constTypesNum_(NULL),
520    nNonLinear_(0),
521    tiny_(source.tiny_),
522    veryTiny_(source.veryTiny_),
523    infty_(source.infty_),
524    exposeWarmStart_(source.exposeWarmStart_),
525    firstSolve_(true),
526    cutStrengthener_(source.cutStrengthener_),
527    oaMessages_(),
528    oaHandler_(NULL),
529    strong_branching_solver_(source.strong_branching_solver_)
530{
531   if(0 && defaultHandler()){
532     messageHandler()->setLogLevel(source.messageHandler()->logLevel());
533   }
534  //Pass in message handler
535  if(0 && source.messageHandler())
536    passInMessageHandler(source.messageHandler());
537   //Copy options from old application
538  if(IsValid(source.tminlp_)) {
539    problem_ = source.problem_->clone();
540    feasibilityProblem_ = new TNLP2FPNLP
541        (SmartPtr<TNLP>(GetRawPtr(problem_)), source.feasibilityProblem_);
542    if(feasibility_mode_)
543      problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
544    else
545      problem_to_optimize_ = GetRawPtr(problem_);
546    pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
547
548    setAuxiliaryInfo(source.getAuxiliaryInfo());
549    // Copy options from old application
550    app_ = source.app_->clone();
551    for(std::list<Ipopt::SmartPtr<TNLPSolver> >::const_iterator i = source.debug_apps_.begin();
552        i != source.debug_apps_.end() ; i++){
553      debug_apps_.push_back((*i)->clone());
554    }
555    testOthers_ = source.testOthers_;
556  }
557  else {
558    throw SimpleError("Don't know how to copy an empty IpoptInterface.",
559        "copy constructor");
560  }
561
562  warmstart_ = source.warmstart_ ? source.warmstart_->clone() : NULL;
563
564  if(source.obj_) {
565    obj_ = new double[source.getNumCols()];
566    CoinCopyN(source.obj_, source.getNumCols(), obj_);
567  }
568
569
570   oaHandler_ = new OaMessageHandler(*source.oaHandler_);;
571
572}
573
574OsiSolverInterface * 
575OsiTMINLPInterface::clone(bool copyData ) const
576{
577  if(copyData)
578    return new OsiTMINLPInterface(*this);
579  else return new OsiTMINLPInterface;
580}
581
582/// Assignment operator
583OsiTMINLPInterface & OsiTMINLPInterface::operator=(const OsiTMINLPInterface& rhs)
584{
585  if(this!= &rhs) {
586    OsiSolverInterface::operator=(rhs);
587    OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
588    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
589    totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
590    totalIterations_ = rhs.totalIterations_;
591    maxRandomRadius_ = rhs.maxRandomRadius_;
592    hasVarNamesFile_ = rhs.hasVarNamesFile_;
593    pushValue_ = rhs.pushValue_;
594
595    delete warmstart_;
596    warmstart_ = NULL;
597
598    if(IsValid(rhs.tminlp_)) {
599
600      tminlp_ = rhs.tminlp_;
601      problem_ = new TMINLP2TNLP(tminlp_);
602      problem_to_optimize_ = GetRawPtr(problem_);
603      app_ = rhs.app_->clone();
604
605      warmstart_ = rhs.warmstart_ ? rhs.warmstart_->clone() : NULL;
606
607      feasibilityProblem_ = new TNLP2FPNLP
608          (SmartPtr<TNLP>(GetRawPtr(problem_)));
609      nnz_jac = rhs.nnz_jac;
610
611      if(constTypes_ != NULL) {
612        delete [] constTypes_;
613        constTypes_ = NULL;
614      }
615      if(rhs.constTypes_ != NULL) {
616        constTypes_ = new TNLP::LinearityType[getNumRows()];
617        CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
618      }
619/*
620      if(constTypesNum_ != NULL) {
621        delete [] constTypesNum_;
622        constTypesNum_ = NULL;
623      }
624      if(rhs.constTypesNum_ != NULL) {
625        constTypesNum_ = new int[getNumRows()];
626        CoinCopyN(rhs.constTypesNum_, getNumRows(), constTypesNum_);
627      }
628*/
629      if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
630        jValues_ = new double [nnz_jac];
631        jCol_    = new Index [nnz_jac];
632        jRow_    = new Index [nnz_jac];
633        CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
634        CoinCopyN(rhs.jCol_    , nnz_jac,jCol_    );
635        CoinCopyN(rhs.jRow_    , nnz_jac,jRow_    );
636      }
637      else if(nnz_jac > 0) {
638        throw CoinError("Arrays for storing jacobian are inconsistant.",
639            "copy constructor",
640            "IpoptOAInterface");
641      }
642      tiny_ = rhs.tiny_;
643      veryTiny_ = rhs.veryTiny_;
644      infty_ = rhs.infty_;
645      exposeWarmStart_ = rhs.exposeWarmStart_;
646
647    }
648    else {
649      tminlp_ =NULL;
650      problem_ = NULL;
651      app_ = NULL;
652      feasibilityProblem_ = NULL;
653    }
654
655
656    if(obj_) {
657      delete [] obj_;
658      obj_ = NULL;
659    }
660    if(rhs.obj_) {
661      obj_ = new double[rhs.getNumCols()];
662      CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
663    }
664
665    hasVarNamesFile_ = rhs.hasVarNamesFile_;
666
667    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
668    totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
669    totalIterations_ = rhs.totalIterations_;
670    maxRandomRadius_ = rhs.maxRandomRadius_;
671    pushValue_ = rhs.pushValue_;
672    numRetryInitial_ = rhs.numRetryInitial_;
673    numRetryResolve_ = rhs.numRetryResolve_;
674    numRetryInfeasibles_ = rhs.numRetryInfeasibles_;
675    numRetryUnsolved_ = rhs.numRetryUnsolved_;
676    pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
677    numIterationSuspect_ = rhs.numIterationSuspect_;
678
679    hasBeenOptimized_ = rhs.hasBeenOptimized_;
680    cutStrengthener_ = rhs.cutStrengthener_;
681
682    delete oaHandler_;
683    oaHandler_ = new OaMessageHandler(*rhs.oaHandler_);
684    strong_branching_solver_ = rhs.strong_branching_solver_;
685
686    freeCachedData();
687  }
688  return *this;
689}
690
691const SmartPtr<OptionsList> OsiTMINLPInterface::options() const
692{
693  if(!IsValid(app_)) {
694    messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
695    return NULL;
696  }
697  else
698    return app_->options();
699}
700
701SmartPtr<OptionsList> OsiTMINLPInterface::options()
702{
703  if(!IsValid(app_)) {
704    messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
705    return NULL;
706  }
707  else
708    return app_->options();
709}
710
711/// Destructor
712OsiTMINLPInterface::~OsiTMINLPInterface ()
713{
714  freeCachedData();
715  delete [] jRow_;
716  delete [] jCol_;
717  delete [] jValues_;
718  delete [] constTypes_;
719  delete [] obj_;
720  delete oaHandler_;
721  delete warmstart_;
722}
723
724void
725OsiTMINLPInterface::freeCachedColRim()
726{
727    if(reducedCosts_!=NULL) {
728    delete []  reducedCosts_;
729    reducedCosts_ = NULL;
730  }
731
732}
733
734void
735OsiTMINLPInterface::freeCachedRowRim()
736{
737  if(rowsense_!=NULL) {
738    delete []  rowsense_;
739    rowsense_ = NULL;
740  }
741  if(rhs_!=NULL) {
742    delete []  rhs_;
743    rhs_ = NULL;
744  }
745  if(rowrange_!=NULL) {
746    delete []  rowrange_;
747    rowrange_ = NULL;
748  }
749  //   if(dualsol_!=NULL)
750  //     {
751  //       delete []  dualsol_; dualsol_ = NULL;
752  //     }
753}
754
755void
756OsiTMINLPInterface::freeCachedData()
757{
758  freeCachedColRim();
759  freeCachedRowRim();
760}
761
762const char * OsiTMINLPInterface::OPT_SYMB="OPT";
763const char * OsiTMINLPInterface::FAILED_SYMB="FAILED";
764const char * OsiTMINLPInterface::UNBOUND_SYMB="UNBOUNDED";
765const char * OsiTMINLPInterface::INFEAS_SYMB="INFEAS";
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,"resolveForCost");
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()<<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,"resolveForRobustness");
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()<<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,"resolveForRobustness");
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()<<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  double * g = new double[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  CoinPackedVector * cuts = new CoinPackedVector[nNonLinear_ + 1];
1766  double * lb = new double[nNonLinear_ + 1];
1767  double * ub = new double[nNonLinear_ + 1];
1768
1769  int * row2cutIdx = new int[m];//store correspondance between index of row and index of cut (some cuts are not generated for rows because linear, or not binding). -1 if constraint does not generate a cut, otherwise index in cuts.
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#if 0
1783      if(fabs(duals[rowIdx]) == 0.)
1784      {
1785        row2cutIdx[rowIdx] = -1;
1786#ifdef NDEBUG
1787        (*handler_)<<"non binding constraint"<<CoinMessageEol;
1788#endif
1789        continue;
1790      }
1791#endif
1792      row2cutIdx[rowIdx] = numCuts;
1793      if(rowLower[rowIdx] > - nlp_infty)
1794        lb[numCuts] = rowLower[rowIdx] - g[rowIdx];
1795      else
1796        lb[numCuts] = - infty;
1797      if(rowUpper[rowIdx] < nlp_infty)
1798        ub[numCuts] = rowUpper[rowIdx] - g[rowIdx];
1799      else
1800        ub[numCuts] = infty;
1801      if(rowLower[rowIdx] > -infty && rowUpper[rowIdx] < infty)
1802      {
1803        if(duals[rowIdx] >= 0)// <= inequality
1804          lb[numCuts] = - infty;
1805        if(duals[rowIdx] <= 0)// >= inequality
1806          ub[numCuts] = infty;
1807      }
1808     
1809      numCuts++;
1810    }
1811    else
1812      row2cutIdx[rowIdx] = -1;
1813  }
1814
1815
1816  for(int i = 0 ; i < nnz_jac_g ; i++) {
1817    const int &rowIdx = jRow_[i];
1818    const int & cutIdx = row2cutIdx[ rowIdx ];
1819    if(cutIdx != -1) {
1820      const int &colIdx = jCol_[i];
1821      //"clean" coefficient
1822      if(cleanNnz(jValues_[i],colLower[colIdx], colUpper[colIdx],
1823                  rowLower[rowIdx], rowUpper[rowIdx],
1824                  x[colIdx],
1825                  lb[cutIdx],
1826                  ub[cutIdx], tiny_, veryTiny_)) {
1827        cuts[cutIdx].insert(colIdx,jValues_[i]);
1828        if(lb[cutIdx] > - infty)
1829          lb[cutIdx] += jValues_[i] * x[colIdx];
1830        if(ub[cutIdx] < infty)
1831          ub[cutIdx] += jValues_[i] * x[colIdx];
1832      }
1833    }
1834  }
1835
1836  int * cut2rowIdx = NULL;
1837  if (IsValid(cutStrengthener_) || oaHandler_->logLevel() > 0) {
1838    cut2rowIdx = new int [numCuts];// Store correspondance between indices of cut and indices of rows. For each cut
1839    for(int rowIdx = 0 ; rowIdx < m ; rowIdx++){
1840       if(row2cutIdx[rowIdx] >= 0){
1841          cut2rowIdx[row2cutIdx[rowIdx]] = rowIdx;
1842       }
1843    }
1844  }
1845
1846  for(int cutIdx = 0; cutIdx < numCuts; cutIdx++) {
1847    //Compute cut violation
1848    if(x2 != NULL) {
1849      double rhs = cuts[cutIdx].dotProduct(x2);
1850      double violation = 0.;
1851      violation = std::max(violation, rhs - ub[cutIdx]);
1852      violation = std::max(violation, lb[cutIdx] - rhs);
1853      if(violation < theta) {
1854        if(oaHandler_->logLevel() > 0)
1855          oaHandler_->message(CUT_NOT_VIOLATED_ENOUGH, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
1856        continue;}
1857        if(oaHandler_->logLevel() > 0)
1858          oaHandler_->message(VIOLATED_OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
1859    }
1860    else if (oaHandler_->logLevel() > 0)
1861      oaHandler_->message(OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<CoinMessageEol;
1862  OsiRowCut newCut;
1863    //    if(lb[i]>-1e20) assert (ub[i]>1e20);
1864
1865    if (IsValid(cutStrengthener_)) {
1866      const int& rowIdx = cut2rowIdx[cutIdx];
1867      bool retval =
1868        cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
1869                                       GetRawPtr(problem_), rowIdx,
1870                                       cuts[cutIdx], lb[cutIdx], ub[cutIdx], g[rowIdx],
1871                                       rowLower[rowIdx], rowUpper[rowIdx],
1872                                       n, x, infty);
1873      if (!retval) {
1874        (*messageHandler()) << "error in cutStrengthener_->ComputeCuts\n";
1875        //exit(-2);
1876      }
1877    }
1878    if(global) {
1879      newCut.setGloballyValidAsInteger(1);
1880    }
1881    newCut.setEffectiveness(99.99e99);
1882    newCut.setLb(lb[cutIdx]);
1883    newCut.setUb(ub[cutIdx]);
1884    newCut.setRow(cuts[cutIdx]);
1885    //    CftValidator validator;
1886    //    validator(newCut);
1887    if(oaHandler_->logLevel()>2){
1888      oaHandler_->print(newCut);}
1889    cs.insert(newCut);
1890  }
1891
1892  delete [] g;
1893  delete [] cuts;
1894  delete [] row2cutIdx;
1895  delete [] cut2rowIdx;
1896
1897  if(getObj == 2 || (getObj && !problem_->hasLinearObjective())) { // Get the objective cuts
1898    double * obj = new double [n];
1899    problem_to_optimize_->eval_grad_f(n, x, 1,obj);
1900    double f;
1901    problem_to_optimize_->eval_f(n, x, 1, f);
1902
1903    CoinPackedVector v;
1904    v.reserve(n);
1905    lb[nNonLinear_] = -f;
1906    ub[nNonLinear_] = -f;
1907    //double minCoeff = 1e50;
1908    for(int i = 0; i<n ; i++) {
1909      if(cleanNnz(obj[i],colLower[i], colUpper[i],
1910          -getInfinity(), 0,
1911          x[i],
1912          lb[nNonLinear_],
1913          ub[nNonLinear_],tiny_, 1e-15)) {
1914        //            minCoeff = min(fabs(obj[i]), minCoeff);
1915        v.insert(i,obj[i]);
1916        lb[nNonLinear_] += obj[i] * x[i];
1917        ub[nNonLinear_] += obj[i] * x[i];
1918      }
1919    }
1920    v.insert(n,-1);
1921    //Compute cut violation
1922    bool genCut = true;
1923    if(x2 != NULL) {
1924      double rhs = v.dotProduct(x2);
1925      double violation = std::max(0., rhs - ub[nNonLinear_]);
1926      if(violation < theta) genCut = false;
1927    }
1928    if(genCut) {
1929      if (IsValid(cutStrengthener_)) {
1930        lb[nNonLinear_] = -infty;
1931        bool retval =
1932          cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
1933                                         GetRawPtr(problem_), -1,
1934                                         v, lb[nNonLinear_], ub[nNonLinear_],
1935                                         ub[nNonLinear_], -infty, 0.,
1936                                         n, x, infty);
1937        if (!retval) {
1938    (*handler_)<< "error in cutStrengthener_->ComputeCuts"<<CoinMessageEol;
1939          //exit(-2);
1940        }
1941      }
1942      OsiRowCut newCut;
1943      if(global)
1944        newCut.setGloballyValidAsInteger(1);
1945      newCut.setEffectiveness(99.99e99);
1946      newCut.setRow(v);
1947      newCut.setLb(-COIN_DBL_MAX/*Infinity*/);
1948      newCut.setUb(ub[nNonLinear_]);
1949      //     CftValidator validator;
1950      //     validator(newCut);
1951      cs.insert(newCut);
1952    }
1953    delete [] obj;
1954    }
1955    else{
1956      printf("Objective is linear\n");
1957    }
1958
1959  delete []lb;
1960  delete[]ub;
1961}
1962
1963
1964
1965/** Get the outer approximation of a single constraint at the point x.
1966*/
1967void
1968OsiTMINLPInterface::getConstraintOuterApproximation(OsiCuts &cs, int rowIdx, 
1969                                                    const double * x, 
1970                                                    const double * x2, bool global)
1971{
1972  double g;
1973  int * indices = new int[getNumCols()];
1974  double * values = new double[getNumCols()];
1975  int nnz;
1976  problem_->eval_grad_gi(getNumCols(), x, 1, rowIdx, nnz, indices, values);
1977  problem_->eval_gi(getNumCols(),x,1, rowIdx, g);
1978
1979  CoinPackedVector cut;
1980  double lb;
1981  double ub;
1982
1983
1984  const double rowLower = getRowLower()[rowIdx];
1985  const double rowUpper = getRowUpper()[rowIdx];
1986  const double * colLower = getColLower();
1987  const double * colUpper = getColUpper();
1988  const double dual = (getRowPrice() + 2 * getNumCols())[rowIdx];
1989  double infty = getInfinity();
1990  double nlp_infty = infty_;
1991 
1992  if(rowLower > - nlp_infty)
1993    lb = rowLower - g;
1994  else
1995    lb = - infty;
1996  if(rowUpper < nlp_infty)
1997    ub = rowUpper - g;
1998  else
1999    ub = infty;
2000  if(rowLower > -infty && rowUpper < infty)
2001  {
2002    if(dual >= 0)// <= inequality
2003      lb = - infty;
2004    if(dual <= 0)// >= inequality
2005      ub = infty;
2006  }
2007
2008  for(int i = 0 ; i < nnz; i++) {
2009     const int &colIdx = indices[i];
2010      //"clean" coefficient
2011      if(cleanNnz(values[i],colLower[colIdx], colUpper[colIdx],
2012                  rowLower, rowUpper,
2013                  x[colIdx],
2014                  lb,
2015                  ub, tiny_, veryTiny_)) {
2016        cut.insert(colIdx,values[i]);
2017        if(lb > - infty)
2018          lb += values[i] * x[colIdx];
2019        if(ub < infty)
2020          ub += values[i] * x[colIdx];
2021    }
2022  }
2023
2024  OsiRowCut newCut;
2025
2026  if(global) {
2027    newCut.setGloballyValidAsInteger(1);
2028  }
2029  newCut.setEffectiveness(99.99e99);
2030  newCut.setLb(lb);
2031  newCut.setUb(ub);
2032  newCut.setRow(cut);
2033  cs.insert(newCut);
2034
2035  delete [] indices;
2036  delete [] values;
2037}
2038
2039void
2040OsiTMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
2041                                            double a, double s, int L){
2042  if(! IsValid(feasibilityProblem_)) {
2043    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2044  }
2045  feasibilityProblem_->set_use_feasibility_pump_objective(true);
2046  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2047  feasibilityProblem_->setLambda(a);
2048  feasibilityProblem_->setSigma(s);
2049  feasibilityProblem_->setNorm(L);
2050  feasibilityProblem_->set_use_cutoff_constraint(false);
2051  feasibilityProblem_->set_use_local_branching_constraint(false); 
2052  problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
2053  feasibility_mode_ = true;
2054}
2055
2056void
2057OsiTMINLPInterface::switchToFeasibilityProblem(int n,const double * x_bar,const int *inds,
2058                                               double rhs_local_branching_constraint){
2059  if(! IsValid(feasibilityProblem_)) {
2060    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2061  }
2062  feasibilityProblem_->set_use_feasibility_pump_objective(false);
2063  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2064  feasibilityProblem_->set_use_cutoff_constraint(false);
2065  feasibilityProblem_->set_use_local_branching_constraint(true); 
2066  feasibilityProblem_->set_rhs_local_branching_constraint(rhs_local_branching_constraint); 
2067  problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
2068  feasibility_mode_ = true;
2069}
2070
2071void
2072OsiTMINLPInterface::switchToOriginalProblem(){
2073  problem_to_optimize_ = GetRawPtr(problem_);
2074  feasibility_mode_ = false;
2075}
2076
2077double
2078OsiTMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
2079                                            double a, double s, int L)
2080{
2081  if(! IsValid(feasibilityProblem_)) {
2082    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2083  }
2084  feasibilityProblem_->set_use_feasibility_pump_objective(true);
2085  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2086  feasibilityProblem_->setLambda(a);
2087  feasibilityProblem_->setSigma(s);
2088  feasibilityProblem_->setNorm(L);
2089  feasibilityProblem_->set_use_cutoff_constraint(false);
2090  feasibilityProblem_->set_use_local_branching_constraint(false); 
2091  nCallOptimizeTNLP_++;
2092  totalNlpSolveTime_-=CoinCpuTime();
2093  SmartPtr<TNLPSolver> app2 = app_->clone();
2094  app2->options()->SetIntegerValue("print_level", (Index) 0);
2095  optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2096  totalNlpSolveTime_+=CoinCpuTime();
2097  hasBeenOptimized_=true;
2098  return getObjValue();
2099}
2100
2101double
2102OsiTMINLPInterface::solveFeasibilityProblem(int n,const double * x_bar,const int *inds, 
2103                                            int L, double cutoff)
2104{
2105  if(! IsValid(feasibilityProblem_)) {
2106    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2107  }
2108  feasibilityProblem_->set_use_feasibility_pump_objective(true);
2109  feasibilityProblem_->set_dist2point_obj(n,(const Number *) x_bar,(const Index *) inds);
2110  feasibilityProblem_->setLambda(1.0);
2111  feasibilityProblem_->setSigma(0.0);
2112  feasibilityProblem_->setNorm(L);
2113  feasibilityProblem_->set_use_cutoff_constraint(true);
2114  feasibilityProblem_->set_cutoff(cutoff);
2115  feasibilityProblem_->set_use_local_branching_constraint(false); 
2116  nCallOptimizeTNLP_++;
2117  totalNlpSolveTime_-=CoinCpuTime();
2118  SmartPtr<TNLPSolver> app2 = app_->clone();
2119  app2->options()->SetIntegerValue("print_level", (Index) 0);
2120  optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2121  totalNlpSolveTime_+=CoinCpuTime();
2122  hasBeenOptimized_=true;
2123  return getObjValue();
2124}
2125
2126#if 0
2127double
2128OsiTMINLPInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs, bool addOnlyViolated, bool global)
2129{
2130  double ret_val = solveFeasibilityProblem(n, x_bar, inds, 1, 0, 2);
2131  getOuterApproximation(cs, getColSolution(), 0, (addOnlyViolated? x_bar:NULL)
2132                        , global);
2133  return ret_val;
2134}
2135#endif
2136
2137static bool WarnedForNonConvexOa=false;
2138
2139void
2140OsiTMINLPInterface::extractLinearRelaxation(OsiSolverInterface &si, 
2141                                            const double * x, bool getObj
2142                                            )
2143{
2144  int n;
2145  int m;
2146  int nnz_jac_g;
2147  int nnz_h_lag;
2148  TNLP::IndexStyleEnum index_style;
2149  //Get problem information
2150  problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2151
2152  //if not allocated allocate spaced for stroring jacobian
2153  //if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2154    initializeJacobianArrays();
2155
2156  //get Jacobian
2157  problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
2158
2159
2160  vector<double> g(m);
2161  problem_to_optimize_->eval_g(n, x, 1, m, g());
2162
2163  vector<double> rowLow(m);
2164  vector<double> rowUp(m);
2165  vector<int> nonBindings(m);;//store non binding constraints (which are to be removed from OA)
2166  int numNonBindings = 0;
2167  const double * rowLower = getRowLower();
2168  const double * rowUpper = getRowUpper();
2169  const double * colLower = getColLower();
2170  const double * colUpper = getColUpper();
2171  const double * duals = getRowPrice() + 2*n;
2172  assert(m==getNumRows());
2173  double infty = si.getInfinity();
2174  double nlp_infty = infty_;
2175 
2176  for(int i = 0 ; i < m ; i++) {
2177    if(constTypes_[i] == TNLP::NON_LINEAR) {
2178      //If constraint is range not binding prepare to remove it
2179      if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty && fabs(duals[i]) == 0.)
2180      {
2181        nonBindings[numNonBindings++] = i;
2182        continue;
2183      }
2184      else
2185        if(rowLower[i] > - nlp_infty){
2186          rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
2187          if(! WarnedForNonConvexOa && rowUpper[i] < nlp_infty){
2188             messageHandler()->message(WARNING_NON_CONVEX_OA, messages_)<<CoinMessageEol;
2189             WarnedForNonConvexOa = true;
2190          }
2191        }
2192      else
2193        rowLow[i] = - infty;
2194      if(rowUpper[i] < nlp_infty)
2195        rowUp[i] =  (rowUpper[i] - g[i]) + 1e-07;
2196      else
2197        rowUp[i] = infty;
2198     
2199      //If equality or ranged constraint only add one side by looking at sign of dual multiplier
2200      if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty)
2201      {
2202        if(duals[i] >= 0.)// <= inequality
2203          rowLow[i] = - infty;
2204        if(duals[i] <= 0.)// >= inequality
2205          rowUp[i] = infty;
2206      }
2207    }
2208    else {
2209      if(rowLower[i] > -nlp_infty){
2210      //   printf("Lower %g ", rowLower[i]);
2211         rowLow[i] = (rowLower[i] - g[i]);
2212      }
2213      else
2214        rowLow[i] = - infty;
2215      if(rowUpper[i] < nlp_infty){
2216      //   printf("Upper %g ", rowUpper[i]);
2217         rowUp[i] =  (rowUpper[i] - g[i]);
2218      }
2219      else
2220        rowUp[i] = infty;
2221    }
2222  }
2223
2224 
2225 
2226  //Then convert everything to a CoinPackedMatrix
2227  //Go through values, clean coefficients and fix bounds
2228  for(int i = 0 ; i < nnz_jac_g ; i++) {
2229    if(constTypes_[jRow_[i]] != TNLP::LINEAR){//For linear just copy is fine.
2230       if(//For other clean tinys
2231       cleanNnz(jValues_[i],colLower[jCol_[i]], colUpper[jCol_[i]],
2232                rowLower[jRow_[i]], rowUpper[jRow_[i]],
2233                x[jCol_[i]],
2234                rowLow[jRow_[i]],
2235                rowUp[jRow_[i]], tiny_, veryTiny_)) {     
2236          rowLow[jRow_[i]] += jValues_[i] * x[jCol_ [i]];
2237          rowUp[jRow_[i]] += jValues_[i] *x[jCol_[i]];
2238       }
2239    }
2240    else {
2241      double value = jValues_[i] * getColSolution()[jCol_[i]];
2242      rowLow[jRow_[i]] += value;
2243      rowUp[jRow_[i]] += value;
2244    } 
2245  }
2246  CoinPackedMatrix mat(true, jRow_, jCol_, jValues_, nnz_jac_g);
2247  mat.setDimensions(m,n); // In case matrix was empty, this should be enough
2248 
2249  //remove non-bindings equality constraints
2250  mat.deleteRows(numNonBindings, nonBindings());
2251
2252  int numcols=getNumCols();
2253  vector<double> obj(numcols);
2254  for(int i = 0 ; i < numcols ; i++)
2255    obj[i] = 0.;
2256 
2257 
2258  si.loadProblem(mat, getColLower(), getColUpper(), obj(), rowLow(), rowUp());
2259  for(int i = 0 ; i < getNumCols() ; i++) {
2260    if(isInteger(i))
2261      si.setInteger(i);
2262  }
2263  if(getObj) {
2264     bool addObjVar = false;
2265     if(problem_->hasLinearObjective()){
2266       double zero;
2267       problem_to_optimize_->eval_f(n, x, 1, zero);
2268       si.setDblParam(OsiObjOffset, -zero);
2269       //Copy the linear objective and don't create a dummy variable.
2270       problem_to_optimize_->eval_grad_f(n, x, 1,obj());
2271       si.setObjective(obj());
2272    }
2273    else {
2274      addObjVar = true;
2275   }
2276   
2277   if(addObjVar){
2278      addObjectiveFunction(si, x);
2279    }
2280  }
2281//  si.writeMpsNative("OA.mps",NULL, NULL, 1);
2282}
2283
2284void 
2285OsiTMINLPInterface::addObjectiveFunction(OsiSolverInterface &si, 
2286                                         const double *x){
2287  const double * colLower = getColLower();
2288  const double * colUpper = getColUpper();
2289  int numcols = getNumCols();
2290  assert(numcols == si.getNumCols() );
2291  vector<double> obj(numcols);
2292      problem_to_optimize_->eval_grad_f(numcols, x, 1,obj());
2293      //add variable alpha
2294      //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
2295      CoinPackedVector a;
2296      si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
2297 
2298      // Now get the objective cuts
2299      // get the gradient, pack it and add the cut
2300      double ub;
2301      problem_to_optimize_->eval_f(numcols, x, 1, ub);
2302      ub*=-1;
2303      double lb = -1e300;
2304      CoinPackedVector objCut;
2305      CoinPackedVector * v = &objCut;
2306      v->reserve(numcols+1);
2307      for(int i = 0; i<numcols ; i++) {
2308       if(si.getNumRows())
2309       {
2310        if(cleanNnz(obj[i],colLower[i], colUpper[i],
2311            -getInfinity(), 0,
2312            x[i],
2313            lb,
2314            ub, tiny_, veryTiny_)) {
2315          v->insert(i,obj[i]);
2316          lb += obj[i] * x[i];
2317          ub += obj[i] * x[i];
2318        }
2319       }
2320       else //Unconstrained problem can not put clean coefficient
2321       {
2322           if(cleanNnz(obj[i],colLower[i], colUpper[i],
2323            -getInfinity(), 0,
2324            x[i],
2325            lb,
2326            ub, 1e-03, 1e-08)) {
2327          v->insert(i,obj[i]);
2328          lb += obj[i] * x[i];
2329          ub += obj[i] * x[i];
2330           }
2331       }
2332      }
2333    v->insert(numcols,-1);
2334    si.addRow(objCut, lb, ub);
2335    }
2336
2337/** Add a collection of linear cuts to problem formulation.*/
2338void 
2339OsiTMINLPInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
2340{ 
2341  if(numberCuts)
2342    freeCachedRowRim();
2343  const OsiRowCut ** cutsPtrs = new const OsiRowCut*[numberCuts];
2344  for(int i = 0 ; i < numberCuts ; i++)
2345  {
2346    cutsPtrs[i] = &cuts[i];
2347  }
2348  problem_->addCuts(numberCuts, cutsPtrs);
2349  delete [] cutsPtrs;
2350}
2351
2352void
2353OsiTMINLPInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
2354    const char * whereFrom)
2355{
2356  totalNlpSolveTime_-=CoinCpuTime();
2357  if(warmStarted)
2358    optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
2359  else
2360    optimizationStatus_ = app_->OptimizeTNLP(GetRawPtr(problem_to_optimize_));
2361  totalNlpSolveTime_+=CoinCpuTime();
2362  nCallOptimizeTNLP_++;
2363  hasBeenOptimized_ = true;
2364 
2365 
2366  //Options should have been printed if not done already turn off Ipopt output
2367  if(!hasPrintedOptions) {
2368    hasPrintedOptions = 1;
2369    //app_->Options()->SetIntegerValue("print_level",0, true, true);
2370    app_->options()->SetStringValue("print_user_options","no", false, true);
2371  }
2372 
2373  bool otherDisagree = false ;
2374#if 0
2375  if(optimizationStatus_ == TNLPSolver::notEnoughFreedom)//Too few degrees of freedom
2376  {
2377    (*messageHandler())<<"Too few degrees of freedom...."<<CoinMessageEol;
2378    int numrows = getNumRows();
2379    int numcols = getNumCols();
2380   
2381    const double * colLower = getColLower();
2382    const double * colUpper = getColUpper();
2383   
2384   
2385    const double * rowLower = getRowLower();
2386    const double * rowUpper = getRowUpper();
2387   
2388    int numberFixed = 0;
2389    for(int i = 0 ; i < numcols ; i++)
2390    {
2391      if(colUpper[i] - colLower[i] <= INT_BIAS)
2392            {
2393              numberFixed++;
2394            }
2395    }
2396    int numberEqualities = 0;
2397    for(int i = 0 ; i < numrows ; i++)
2398    {
2399      if(rowUpper[i] - rowLower[i] <= INT_BIAS)
2400            {
2401              numberEqualities++;
2402            }     
2403    }
2404    if(numcols - numberFixed > numberEqualities || numcols < numberEqualities)
2405    {
2406      std::string probName;
2407      getStrParam(OsiProbName, probName);
2408      throw newUnsolvedError(app_->errorCode(), problem_, probName);
2409    }
2410    double * saveColLow = CoinCopyOfArray(getColLower(), getNumCols());
2411    double * saveColUp = CoinCopyOfArray(getColUpper(), getNumCols());
2412   
2413    for(int i = 0; i < numcols && numcols - numberFixed <= numberEqualities ; i++)
2414    {
2415      if(colUpper[i] - colLower[i] <= INT_BIAS)
2416            {
2417              setColLower(i, saveColLow[i]-1e-06);
2418              setColUpper(i, saveColLow[i]+1e-06);
2419              numberFixed--;
2420            }
2421    }
2422    solveAndCheckErrors(warmStarted, throwOnFailure, whereFrom);
2423    //restore
2424    for(int i = 0; i < numcols && numcols - numberFixed < numrows ; i++)
2425    {
2426      problem_->SetVariableLowerBound(i,saveColLow[i]);
2427      problem_->SetVariableUpperBound(i,saveColUp[i]);
2428    }
2429    delete [] saveColLow;
2430    delete [] saveColUp;
2431    return;
2432  }
2433  else
2434#endif
2435    if(!app_->isRecoverable(optimizationStatus_))//Solver failed and the error can not be recovered, throw it
2436    {
2437      std::string probName;
2438      getStrParam(OsiProbName, probName);
2439      throw newUnsolvedError(app_->errorCode(), problem_, probName);
2440    }
2441    else if(testOthers_ && !app_->isError(optimizationStatus_)){
2442      Ipopt::SmartPtr<TMINLP2TNLP> problem_copy = problem_->clone();
2443      //Try other solvers and see if they agree
2444      int f =1;
2445      for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2446          i != debug_apps_.end() ; i++){
2447        TNLPSolver::ReturnStatus otherStatus = (*i)->OptimizeTNLP(GetRawPtr(problem_copy));
2448       messageHandler()->message(LOG_LINE, messages_)
2449         <<'d'<<f++<<statusAsString(otherStatus)<<problem_copy->obj_value()
2450         <<(*i)->IterationCount()<<(*i)->CPUTime()<<CoinMessageEol;
2451        if(!(*i)->isError(otherStatus)){
2452           CoinRelFltEq eq(1e-05);
2453           if(otherStatus != optimizationStatus_){
2454             otherDisagree = true;
2455             messageHandler()->message(SOLVER_DISAGREE_STATUS, messages_)
2456             <<app_->solverName()<<statusAsString()
2457             <<(*i)->solverName()<<statusAsString(otherStatus)<<CoinMessageEol; 
2458           }
2459           else if(isProvenOptimal() && !eq(problem_->obj_value(),problem_copy->obj_value()))
2460           {
2461             otherDisagree = true;
2462             messageHandler()->message(SOLVER_DISAGREE_VALUE, messages_)
2463             <<app_->solverName()<<problem_->obj_value()
2464             <<(*i)->solverName()<<problem_copy->obj_value()<<CoinMessageEol; 
2465           }
2466        }
2467     }
2468  }
2469  else if(app_->isError(optimizationStatus_) && ! debug_apps_.empty()){
2470      int f =1;
2471      for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2472          i != debug_apps_.end() && app_->isError(optimizationStatus_) ; i++){
2473        optimizationStatus_ = (*i)->OptimizeTNLP(GetRawPtr(problem_));
2474        messageHandler()->message(LOG_LINE, messages_)
2475          <<'d'<<f++<<statusAsString(optimizationStatus_)<<problem_->obj_value()
2476          <<(*i)->IterationCount()<<(*i)->CPUTime()<<CoinMessageEol;
2477      }
2478  }
2479  try{
2480    totalIterations_ += app_->IterationCount();
2481  }
2482  catch(SimpleError &E)
2483  {
2484    if (throwOnFailure)//something failed throw
2485    {
2486      throw SimpleError("No statistics available from Ipopt",whereFrom);
2487    }
2488    else {
2489      return;
2490    }
2491  }
2492  if(problem_->hasUpperBoundingObjective()){//Check if solution is integer and recompute objective value using alternative objective function
2493    const double * sol = getColSolution();
2494    bool integerSol = true;
2495    double intTol = 1e-08;
2496    if(objects()){
2497      int nObjects = numberObjects();
2498      OsiObject ** object = objects();
2499      for(int i = 0 ; i< nObjects ; i++){
2500        int dummy;
2501        if(object[i]->infeasibility(this,dummy) > intTol)
2502        {
2503          integerSol=false;
2504          break;
2505        }
2506      }
2507    }
2508    else{//Only works for integer constraints
2509      int numcols = getNumCols();
2510      for(int i = 0 ; i < numcols ; i++){
2511        if(isInteger(i) || isBinary(i)){
2512          if(fabs(sol[i] - floor(sol[i]+0.5)) > intTol){
2513            integerSol = false;
2514            break;
2515          }
2516        }
2517      }
2518    }
2519    if(integerSol&&isProvenOptimal()){
2520      double help= problem_->evaluateUpperBoundingFunction(sol);
2521     
2522
2523      OsiAuxInfo * auxInfo = getAuxiliaryInfo();
2524      Bonmin::AuxInfo * bonInfo = dynamic_cast<Bonmin::AuxInfo *>(auxInfo);
2525      if(bonInfo!=0)
2526      {
2527       
2528        if(help<bonInfo->bestObj2())
2529        {
2530          bonInfo->setBestObj2(help);
2531          bonInfo->setBestSolution2(getNumCols(), const_cast<double *>(getColSolution()));
2532
2533           messageHandler()->message(ALTERNATE_OBJECTIVE, messages_)
2534           <<help<<CoinMessageEol;
2535        }
2536      }
2537      else {
2538        printf("\nWARNING: the algorithm selected does not consider the second objective function\n");
2539      }
2540    }
2541  }
2542  messageHandler()->message(IPOPT_SUMMARY, messages_)
2543    <<whereFrom<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
2544 
2545  if((nCallOptimizeTNLP_ % 20) == 1)
2546    messageHandler()->message(LOG_HEAD, messages_)<<CoinMessageEol;
2547 
2548 
2549  if ( (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) ||
2550       ( otherDisagree )){
2551    messageHandler()->message(SUSPECT_PROBLEM,
2552                              messages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
2553    std::string subProbName;
2554    getStrParam(OsiProbName, subProbName);
2555    std::ostringstream os;
2556    os<<"_"<<nCallOptimizeTNLP_;
2557    subProbName+=os.str();
2558    problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
2559  }
2560 
2561}
2562
2563////////////////////////////////////////////////////////////////////
2564// Solve Methods                                                  //
2565////////////////////////////////////////////////////////////////////
2566/// Solve initial continuous relaxation
2567void OsiTMINLPInterface::initialSolve()
2568{
2569  assert(IsValid(app_));
2570  assert(IsValid(problem_));
2571
2572  // Discard warmstart_ if we had one
2573  delete warmstart_;
2574  warmstart_ = NULL;
2575 
2576  if(!hasPrintedOptions) {
2577    int printOptions;
2578    app_->options()->GetEnumValue("print_user_options",printOptions,app_->prefix());
2579    if(printOptions)
2580      app_->options()->SetStringValue("print_user_options","yes",true,true);
2581  }
2582  if(exposeWarmStart_)
2583    app_->disableWarmStart(); 
2584  solveAndCheckErrors(0,1,"initialSolve");
2585 
2586  //Options should have been printed if not done already turn off Ipopt output
2587  if(!hasPrintedOptions) {
2588    hasPrintedOptions = 1;
2589    app_->options()->SetStringValue("print_user_options","no");
2590    app_->options()->SetIntegerValue("print_level",0);
2591  }
2592 
2593  messageHandler()->message(LOG_FIRST_LINE, messages_)<<nCallOptimizeTNLP_
2594                                                      <<statusAsString()
2595                                                      <<getObjValue()
2596                                                      <<app_->IterationCount()
2597                                                      <<app_->CPUTime()
2598                                                      <<CoinMessageEol;
2599 
2600  int numRetry = firstSolve_ ? numRetryInitial_ : numRetryResolve_;
2601  if(isAbandoned()) {
2602    resolveForRobustness(numRetryUnsolved_);
2603  }
2604  else if(numRetry)
2605    {
2606      resolveForCost(numRetry, numRetryInitial_ > 0);
2607      /** Only do it once at the root.*/
2608      numRetryInitial_ = 0;
2609    }
2610  firstSolve_ = false;
2611
2612  // if warmstart_ is not empty then had to use resolveFor... and that created
2613  // the warmstart at the end, and we have nothing to do here. Otherwise...
2614  if (! warmstart_ && ! isAbandoned()) {
2615    if (exposeWarmStart_) {
2616      warmstart_ = app_->getWarmStart(problem_);
2617    }
2618  }
2619}
2620
2621/** Resolve the continuous relaxation after problem modification.
2622 * \note for Ipopt, same as resolve */
2623void
2624OsiTMINLPInterface::resolve()
2625{
2626  assert(IsValid(app_));
2627  assert(IsValid(problem_));
2628 
2629  int has_warmstart = warmstart_ == NULL ? 0 : 1;
2630  if(warmstart_ == NULL) has_warmstart = 0;
2631  else if(!app_->warmStartIsValid(warmstart_)) has_warmstart = 1;
2632  else has_warmstart = 2;
2633  if (has_warmstart < 2) {
2634    // empty (or unrecognized) warmstart
2635    initialSolve();
2636    return;
2637  }
2638  app_->setWarmStart(warmstart_, problem_);
2639  delete warmstart_;
2640  warmstart_ = NULL;
2641
2642  if (INT_BIAS > 0.) {
2643    app_->options()->SetStringValue("warm_start_same_structure", "yes");
2644  }
2645  else {
2646    app_->options()->SetStringValue("warm_start_same_structure", "no");
2647  }
2648
2649  if(problem_->duals_init() != NULL)
2650    app_->enableWarmStart();
2651  else app_->disableWarmStart();
2652  solveAndCheckErrors(1,1,"resolve");
2653 
2654  messageHandler()->message(LOG_FIRST_LINE, messages_)<<nCallOptimizeTNLP_
2655                                                      <<statusAsString()
2656                                                      <<getObjValue()
2657                                                      <<app_->IterationCount()
2658                                                      <<app_->CPUTime()<<CoinMessageEol;
2659 
2660  if(isAbandoned()) {
2661    resolveForRobustness(numRetryUnsolved_);
2662  }
2663  else if(numRetryResolve_ ||
2664          (numRetryInfeasibles_ && isProvenPrimalInfeasible() ))
2665    resolveForCost(std::max(numRetryResolve_, numRetryInfeasibles_), 0);
2666
2667  // if warmstart_ is not empty then had to use resolveFor... and that created
2668  // the warmstart at the end, and we have nothing to do here. Otherwise...
2669  if (! warmstart_ && ! isAbandoned()) {
2670    if (exposeWarmStart_) {
2671      warmstart_ = app_->getWarmStart(problem_);
2672    }
2673  }
2674}
2675
2676
2677////////////////////////////////////////////////////////////////
2678// Methods returning info on how the solution process terminated  //
2679////////////////////////////////////////////////////////////////
2680/// Are there a numerical difficulties?
2681bool OsiTMINLPInterface::isAbandoned() const
2682{
2683  return (
2684        (optimizationStatus_==TNLPSolver::iterationLimit)||
2685        (optimizationStatus_==TNLPSolver::computationError)||
2686        (optimizationStatus_==TNLPSolver::illDefinedProblem)||
2687        (optimizationStatus_==TNLPSolver::illegalOption)||
2688        (optimizationStatus_==TNLPSolver::externalException)|
2689        (optimizationStatus_==TNLPSolver::exception)
2690      );
2691}
2692
2693/// Is optimality proven?
2694bool OsiTMINLPInterface::isProvenOptimal() const
2695{
2696  return (optimizationStatus_==TNLPSolver::solvedOptimal) ||
2697          (optimizationStatus_==TNLPSolver::solvedOptimalTol);
2698}
2699/// Is primal infeasiblity proven?
2700bool OsiTMINLPInterface::isProvenPrimalInfeasible() const
2701{
2702  return (optimizationStatus_ == TNLPSolver::provenInfeasible);
2703}
2704/// Is dual infeasiblity proven?
2705bool OsiTMINLPInterface::isProvenDualInfeasible() const
2706{
2707  return (optimizationStatus_ == TNLPSolver::unbounded);
2708}
2709/// Is the given primal objective limit reached?
2710bool OsiTMINLPInterface::isPrimalObjectiveLimitReached() const
2711{
2712  (*handler_)<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2713  return 0;
2714}
2715/// Is the given dual objective limit reached?
2716bool OsiTMINLPInterface::isDualObjectiveLimitReached() const
2717{
2718  //  (*messageHandler_)<<"Warning : isDualObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2719  return (optimizationStatus_==TNLPSolver::unbounded);
2720
2721}
2722/// Iteration limit reached?
2723bool OsiTMINLPInterface::isIterationLimitReached() const
2724{
2725  return (optimizationStatus_==TNLPSolver::iterationLimit);
2726}
2727
2728void
2729OsiTMINLPInterface::extractInterfaceParams()
2730{
2731  if (IsValid(app_)) {
2732    int logLevel;
2733    app_->options()->GetIntegerValue("nlp_log_level", logLevel,app_->prefix());
2734    messageHandler()->setLogLevel(logLevel);
2735
2736#ifdef COIN_HAS_FILTERSQP
2737    FilterSolver * filter = dynamic_cast<FilterSolver *>(GetRawPtr(app_));
2738
2739    bool is_given =
2740#endif
2741      app_->options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,app_->prefix());
2742
2743#ifdef COIN_HAS_FILTERSQP
2744    if(filter && !is_given){
2745      // overwriting default for filter
2746      maxRandomRadius_ = 10.;
2747    }
2748#endif
2749   
2750   int oaCgLogLevel = 0;
2751   app_->options()->GetIntegerValue("oa_cuts_log_level", oaCgLogLevel,app_->prefix());
2752   oaHandler_->setLogLevel(oaCgLogLevel); 
2753   
2754    int exposeWs = false;
2755    app_->options()->GetEnumValue("warm_start", exposeWs, app_->prefix());
2756    setExposeWarmStart(exposeWs > 0);
2757   
2758    app_->options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,app_->prefix());
2759    app_->options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,app_->prefix());
2760    app_->options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,app_->prefix());
2761    app_->options()->GetIntegerValue("num_resolve_at_infeasibles", numRetryInfeasibles_,app_->prefix());
2762    app_->options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,app_->prefix());
2763    app_->options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,app_->prefix());
2764    app_->options()->GetNumericValue
2765    ("warm_start_bound_frac" ,pushValue_,app_->prefix());
2766    app_->options()->GetNumericValue("tiny_element",tiny_,app_->prefix());
2767    app_->options()->GetNumericValue("very_tiny_element",veryTiny_,app_->prefix());
2768    app_->options()->GetNumericValue("random_point_perturbation_interval",max_perturbation_,app_->prefix());
2769    app_->options()->GetEnumValue("random_point_type",randomGenerationType_,app_->prefix());
2770    int cut_strengthening_type;
2771    app_->options()->GetEnumValue("cut_strengthening_type", cut_strengthening_type,app_->prefix());
2772
2773    if (cut_strengthening_type != CS_None) {
2774      // TNLP solver to be used in the cut strengthener
2775      cutStrengthener_ = new CutStrengthener(app_->clone(), app_->options());
2776    }
2777  }
2778}
2779
2780void
2781OsiTMINLPInterface::SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)
2782{
2783  strong_branching_solver_ = strong_branching_solver;
2784}
2785
2786  //#define STRONG_COMPARE
2787#ifdef STRONG_COMPARE
2788  static double objorig;
2789#endif
2790
2791void
2792OsiTMINLPInterface::markHotStart()
2793{
2794  if (IsValid(strong_branching_solver_)) {
2795#ifdef STRONG_COMPARE
2796    // AWDEBUG
2797    OsiSolverInterface::markHotStart();
2798    objorig = getObjValue();
2799#endif
2800    optimizationStatusBeforeHotStart_ = optimizationStatus_;
2801    strong_branching_solver_->markHotStart(this);
2802  }
2803  else {
2804    // Default Implementation
2805    OsiSolverInterface::markHotStart();
2806  }
2807}
2808
2809void
2810OsiTMINLPInterface::solveFromHotStart()
2811{
2812  if (IsValid(strong_branching_solver_)) {
2813#ifdef STRONG_COMPARE
2814    // AWDEBUG
2815    OsiSolverInterface::solveFromHotStart();
2816    double obj_nlp = getObjValue() - objorig;
2817#endif
2818    optimizationStatus_ = strong_branching_solver_->solveFromHotStart(this);
2819    hasBeenOptimized_ = true;
2820#ifdef STRONG_COMPARE
2821    double obj_other = getObjValue() - objorig;
2822    printf("AWDEBUG: Strong Branching results: NLP = %15.8e Other = %15.8e\n",
2823           obj_nlp, obj_other);
2824#endif
2825  }
2826  else {
2827    // Default Implementation
2828    OsiSolverInterface::solveFromHotStart();
2829  }
2830}
2831
2832void
2833OsiTMINLPInterface::unmarkHotStart()
2834{
2835  if (IsValid(strong_branching_solver_)) {
2836#ifdef STRONG_COMPARE
2837    // AWDEBUG
2838    OsiSolverInterface::unmarkHotStart();
2839#endif
2840    strong_branching_solver_->unmarkHotStart(this);
2841    optimizationStatus_ = optimizationStatusBeforeHotStart_;
2842  }
2843  else {
2844    // Default Implementation
2845    OsiSolverInterface::unmarkHotStart();
2846  }
2847}
2848
2849const double * OsiTMINLPInterface::getObjCoefficients() const
2850{
2851  const int n = getNumCols();
2852  delete [] obj_;
2853  obj_ = NULL;
2854  obj_ = new double[n];
2855
2856  bool new_x = true;
2857  const double* x_sol = problem_->x_sol();
2858  bool retval = problem_->eval_grad_f(n, x_sol, new_x, obj_);
2859 
2860  if (!retval) {
2861    // Let's see if that happens - it will cause a crash
2862    fprintf(stderr, "ERROR WHILE EVALUATING GRAD_F in OsiTMINLPInterface::getObjCoefficients()\n");
2863    delete [] obj_;
2864    obj_ = NULL;
2865  }
2866
2867  return obj_;
2868}
2869
2870  /** Sets the TMINLP2TNLP to be used by the interface.*/
2871  void 
2872  OsiTMINLPInterface::use(Ipopt::SmartPtr<TMINLP2TNLP> tminlp2tnlp){
2873     problem_ = tminlp2tnlp;
2874     feasibilityProblem_->use(GetRawPtr(tminlp2tnlp));}
2875
2876}/** end namespace Bonmin*/
2877
Note: See TracBrowser for help on using the repository browser.