source: stable/1.1/Bonmin/src/Interfaces/BonOsiTMINLPInterface.cpp @ 1521

Last change on this file since 1521 was 1521, checked in by pbonami, 10 years ago

Try to fix problem with CPLEX bounds

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