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

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

Remove debug output

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