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

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

Port Claudia's code to trunk

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