source: branches/split/Bonmin/src/Interfaces/BonOsiTMINLPInterface.cpp @ 1739

Last change on this file since 1739 was 1739, checked in by stefan, 9 years ago

COIN_HAS_CLP and COIN_HAS_IPOPT do not exist anymore, and cannot build Bonmin without them anyway

File size: 94.5 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",127);
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",127);
67
68  roptions->AddStringOption2("file_solution",
69       "Write a file bonmin.sol with the solution",
70       "no",
71       "yes","",
72       "no","","");
73  roptions->setOptionExtraInfo("file_solution",127);
74
75  roptions->AddStringOption3("warm_start",
76      "Select the warm start method",
77      "none",
78      "none","No warm start",
79      "optimum","Warm start with direct parent optimum",
80      "interior_point","Warm start with an interior point of direct parent",
81      "This will affect the function getWarmStart(), and as a consequence the warm starting in the various algorithms.");
82  roptions->setOptionExtraInfo("warm_start",8);
83
84  roptions->SetRegisteringCategory("Nlp solution robustness", RegisteredOptions::BonminCategory);
85
86  roptions->AddLowerBoundedNumberOption("max_random_point_radius",
87      "Set max value r for coordinate of a random point.",
88      0.,1,1e5,
89      "When picking a random point, coordinate i will be in the interval [min(max(l,-r),u-r), max(min(u,r),l+r)] "
90      "(where l is the lower bound for the variable and u is its upper bound)");
91  roptions->setOptionExtraInfo("max_random_point_radius",8);
92
93  roptions->AddStringOption3("random_point_type","method to choose a random starting point",
94                             "Jon",
95                             "Jon", "Choose random point uniformly between the bounds",
96                             "Andreas", "perturb the starting point of the problem within a prescribed interval",
97                             "Claudia", "perturb the starting point using the perturbation radius suffix information",
98                             "");
99  roptions->setOptionExtraInfo("random_point_type",8);
100
101    roptions->AddLowerBoundedNumberOption("random_point_perturbation_interval",
102                                           "Amount by which starting point is perturbed when choosing to pick random point by perturbating starting point",
103                                           0.,true, 1.,
104                                           "");
105  roptions->setOptionExtraInfo("random_point_perturbation_interval",8);
106                                           
107
108  roptions->AddLowerBoundedIntegerOption
109  ("num_iterations_suspect",
110   "Number of iterations over which a node is considered \"suspect\" (for debugging purposes only, see detailed documentation).",
111   -1,-1,
112   "When the number of iterations to solve a node is above this number, the subproblem at this"
113   " node is considered to be suspect and it will be outputed in a file (set to -1 to deactivate this).");
114  roptions->setOptionExtraInfo("num_iterations_suspect",127);
115
116 
117
118  roptions->AddLowerBoundedIntegerOption("num_retry_unsolved_random_point",
119      "Number $k$ of times that the algorithm will try to resolve an unsolved NLP with a random starting point "
120      "(we call unsolved an NLP for which Ipopt is not "
121      "able to guarantee optimality within the specified tolerances).",
122      0,0,
123      "When Ipopt fails to solve a continuous NLP sub-problem, if $k > 0$, the algorithm will "
124      "try again to solve the failed NLP with $k$ new randomly chosen starting points "
125      " or until the problem is solved with success.");
126  roptions->setOptionExtraInfo("num_retry_unsolved_random_point",127);
127
128
129  roptions->SetRegisteringCategory("Options for non-convex problems", RegisteredOptions::BonminCategory);
130
131
132  roptions->AddLowerBoundedIntegerOption("num_resolve_at_root",
133      "Number $k$ of tries to resolve the root node with different starting points.",
134      0,0,
135      "The algorithm will solve the root node with $k$ random starting points"
136      " and will keep the best local optimum found.");
137  roptions->setOptionExtraInfo("num_resolve_at_root",8);
138
139  roptions->AddLowerBoundedIntegerOption("num_resolve_at_node",
140      "Number $k$ of tries to resolve a node (other than the root) of the tree with different starting point.",
141      0,0,
142      "The algorithm will solve all the nodes with $k$ different random starting points "
143      "and will keep the best local optimum found.");
144  roptions->setOptionExtraInfo("num_resolve_at_node",8);
145
146  roptions->AddLowerBoundedIntegerOption("num_resolve_at_infeasibles",
147      "Number $k$ of tries to resolve an infeasible node (other than the root) of the tree with different starting point.",
148      0,0,
149      "The algorithm will solve all the infeasible nodes with $k$ different random starting points "
150      "and will keep the best local optimum found.");
151  roptions->setOptionExtraInfo("num_resolve_at_infeasibles",8);
152
153
154  roptions->AddStringOption2("dynamic_def_cutoff_decr",
155      "Do you want to define the parameter cutoff_decr dynamically?",
156      "no",
157      "no", "No, define it statically",
158      "yes","Yes, define it dynamically");
159  roptions->setOptionExtraInfo("dynamic_def_cutoff_decr",8);
160
161  roptions->AddLowerBoundedNumberOption("coeff_var_threshold",
162      "Coefficient of variation threshold (for dynamic definition of cutoff_decr).",
163      0.0,
164      false,
165      0.1);
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  roptions->setOptionExtraInfo("first_perc_for_cutoff_decr",8);
172
173  roptions->AddNumberOption("second_perc_for_cutoff_decr",
174      "The percentage used when, the coeff of variance is greater than the threshold, to compute the cutoff_decr dynamically.",
175      -0.05);
176  roptions->setOptionExtraInfo("second_perc_for_cutoff_decr",8);
177
178  }
179
180static void register_OA_options
181(SmartPtr<RegisteredOptions> roptions)
182{
183        //
184
185  roptions->SetRegisteringCategory("Outer Approximations strengthening", RegisteredOptions::UndocumentedCategory);
186  roptions->AddStringOption2("disjunctive_cut_type",
187      "Determine if and what kind of disjunctive cuts should be computed.",
188      "none",
189      "none", "No disjunctive cuts.",
190      "most-fractional", "If discrete variables present, compute disjunction for most-fractional variable");
191  roptions->setOptionExtraInfo("disjunctive_cut_type",119);
192
193
194 
195 
196  roptions->AddStringOption4("cut_strengthening_type",
197                             "Determines if and what kind of cut strengthening should be performed.",
198                             "none",
199                             "none", "No strengthening of cuts.",
200                             "sglobal", "Strengthen global cuts.",
201                             "uglobal-slocal", "Unstrengthened global and strengthened local cuts",
202                             "sglobal-slocal", "Strengthened global and strengthened local cuts",
203                             "");
204  roptions->setOptionExtraInfo("cut_strengthening_type",119);
205
206  roptions->SetRegisteringCategory("Outer Approximation cuts generation", RegisteredOptions::BonminCategory);
207
208  roptions->AddStringOption2("oa_cuts_scope","Specify if OA cuts added are to be set globally or locally valid",
209                             "global",
210                             "local","Cuts are treated as locally valid",
211                             "global", "Cuts are treated as globally valid",
212                             "");
213  roptions->setOptionExtraInfo("oa_cuts_scope",119);
214
215  roptions->AddStringOption2("add_only_violated_oa","Do we add all OA cuts or only the ones violated by current point?",
216                             "no",
217                             "no","Add all cuts",
218                             "yes","Add only violated Cuts","");
219  roptions->setOptionExtraInfo("add_only_violated_oa",119);
220
221 
222  roptions->AddLowerBoundedNumberOption("tiny_element","Value for tiny element in OA cut",
223      -0.,0,1e-08,
224      "We will remove \"cleanly\" (by relaxing cut) an element lower"
225      " than this.");
226  roptions->setOptionExtraInfo("tiny_element",119);
227
228  roptions->AddLowerBoundedNumberOption("very_tiny_element","Value for very tiny element in OA cut",
229      -0.,0,1e-17,
230      "Algorithm will take the risk of neglecting an element lower"
231      " than this.");
232  roptions->setOptionExtraInfo("very_tiny_element",119);
233
234  roptions->AddLowerBoundedIntegerOption("oa_cuts_log_level",
235                                         "level of log when generating OA cuts.",
236                                         0, 0,
237                                         "0: outputs nothing,\n"
238                                         "1: when a cut is generated, its violation and index of row from which it originates,\n"
239                                         "2: always output violation of the cut.\n"
240                                         "3: output generated cuts incidence vectors.");
241  roptions->setOptionExtraInfo("oa_cuts_log_level",119);
242
243}
244
245
246///Register options
247void
248OsiTMINLPInterface::registerOptions
249(SmartPtr<RegisteredOptions> roptions)
250{
251  // We try to register the options - if those have been registered
252  // already, we catch the exception and don't need to do it again
253  try {
254    register_general_options(roptions);
255    register_OA_options(roptions);
256#ifdef COIN_HAS_FILTERSQP
257    FilterSolver::RegisterOptions(roptions);
258#endif
259    IpoptSolver::RegisterOptions(roptions);
260  }   
261  catch(RegisteredOptions::OPTION_ALREADY_REGISTERED) {
262    // skipping
263  }
264
265
266}
267
268/** To set some application specific defaults. */
269void 
270OsiTMINLPInterface::setAppDefaultOptions(SmartPtr<OptionsList> Options)
271{}
272
273OsiTMINLPInterface::Messages::Messages
274():CoinMessages((int)OSITMINLPINTERFACE_DUMMY_END)
275{
276  strcpy(source_ ,"NLP");
277  ADD_MSG(SOLUTION_FOUND, std_m, 2,
278          "After %d tries found a solution of %g (previous best %g).");
279
280  ADD_MSG(INFEASIBLE_SOLUTION_FOUND, std_m, 2,
281          "After %d tries found an solution of %g infeasible problem.");
282
283  ADD_MSG(UNSOLVED_PROBLEM_FOUND, std_m, 2,
284          "After %d tries found an solution of %g unsolved problem.");
285  ADD_MSG(WARN_SUCCESS_WS, warn_m, 2,
286       "Problem not solved with warm start but solved without");
287
288  ADD_MSG(WARNING_RESOLVING, warn_m,2,
289       "Trying to resolve NLP with different starting point (%d attempts).");
290  ADD_MSG(WARN_SUCCESS_RANDOM, warn_m, 1,
291       "Problem initially not solved but solved with a random starting point (success on %d attempt)");
292  ADD_MSG(WARN_CONTINUING_ON_FAILURE, warn_m, 1,
293       "Warning : continuing branching, while there are unrecovered failures at nodes");
294
295  ADD_MSG(SUSPECT_PROBLEM, std_m, 2, "NLP number %d is suspect (see bounds and start file)");
296  ADD_MSG(IPOPT_SUMMARY, std_m, 2, "Ipopt return (for %s): status %2d, iter count %4d, time %g");
297  ADD_MSG(BETTER_SOL, std_m, 2, "Solution of value %g found on %d'th attempt");
298
299  ADD_MSG(LOG_HEAD, std_m, 1,
300          "\n          "
301          "    Num      Status      Obj             It       time                 Location");
302  ADD_MSG(LOG_LINE, std_m, 1, 
303          "%c    %-8d %-11s %-15.5g %-8d %-8g %20s");
304
305  ADD_MSG(ALTERNATE_OBJECTIVE, std_m, 1, "Objective value recomputed with alternate objective: %g.");
306 
307  ADD_MSG(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, warn_m, 1, 
308         "resolve called before any call to initialSol  can not use warm starts.");
309  ADD_MSG(ERROR_NO_TNLPSOLVER, warn_m, 1,"Can not parse options when no IpApplication has been created");
310  ADD_MSG(WARNING_NON_CONVEX_OA, warn_m, 1,
311          "OA on non-convex constraint is very experimental.");                         
312  ADD_MSG(SOLVER_DISAGREE_STATUS, warn_m, 1, "%s says problem %s, %s says %s.");
313  ADD_MSG(SOLVER_DISAGREE_VALUE, warn_m, 1, "%s gives objective %.16g, %s gives %.16g.");
314
315}
316
317
318void 
319OsiTMINLPInterface::OaMessageHandler::print(OsiRowCut &row){
320  FILE * fp = filePointer();
321  const int &n = row.row().getNumElements();
322  fprintf(fp,"Row cut has %d elements. Lower bound: %g, upper bound %g.\n", n, row.lb(), row.ub());
323  const int * idx = row.row().getIndices();
324  const double * val = row.row().getElements();
325  for(int i = 0 ; i < n ; i++){
326    fprintf(fp,"%g, x%d",val[i], idx[i]);
327    if(i && i % 7 == 0)
328      fprintf(fp,"\n");
329  } 
330}
331
332OsiTMINLPInterface::OaMessages::OaMessages(): CoinMessages((int) OA_MESSAGES_DUMMY_END){
333   strcpy(source_,"OaCg");
334   ADD_MSG(VIOLATED_OA_CUT_GENERATED, std_m, 1,"Row %d, cut violation is %g: Outer approximation cut generated.");
335
336   ADD_MSG(CUT_NOT_VIOLATED_ENOUGH, std_m, 2,"Row %d, cut violation is %g: Outer approximation cut not generated.");
337
338   ADD_MSG(OA_CUT_GENERATED, std_m, 1,"Row %d: Outer approximation cut not generated.");
339}
340bool OsiTMINLPInterface::hasPrintedOptions=0;
341
342////////////////////////////////////////////////////////////////////
343// Constructors and desctructors                                  //
344////////////////////////////////////////////////////////////////////
345/// Default Constructor
346OsiTMINLPInterface::OsiTMINLPInterface():
347    OsiSolverInterface(),
348    tminlp_(NULL),
349    problem_(NULL),
350    problem_to_optimize_(NULL),
351    feasibility_mode_(false),
352    app_(NULL),
353    debug_apps_(),
354    testOthers_(false),
355    warmstart_(NULL),
356    rowsense_(NULL),
357    rhs_(NULL),
358    rowrange_(NULL),
359    reducedCosts_(NULL),
360    OsiDualObjectiveLimit_(1e200),
361    hasVarNamesFile_(true),
362    nCallOptimizeTNLP_(0),
363    totalNlpSolveTime_(0),
364    totalIterations_(0),
365    maxRandomRadius_(1e08),
366    randomGenerationType_(0),
367    max_perturbation_(COIN_DBL_MAX),
368    pushValue_(1e-02),
369    numRetryInitial_(-1),
370    numRetryResolve_(-1),
371    numRetryInfeasibles_(-1),
372    numRetryUnsolved_(1),
373    dynamicCutOff_(0),
374    coeff_var_threshold_(0.1),
375    first_perc_for_cutoff_decr_(-0.02),
376    second_perc_for_cutoff_decr_(-0.05),
377    messages_(),
378    pretendFailIsInfeasible_(0),
379    pretendSucceededNext_(0),
380    hasContinuedAfterNlpFailure_(false),
381    numIterationSuspect_(-1),
382    hasBeenOptimized_(false),
383    obj_(NULL),
384    feasibilityProblem_(NULL),
385    jRow_(NULL),
386    jCol_(NULL),
387    jValues_(NULL),
388    nnz_jac(0),
389    constTypes_(NULL),
390    nNonLinear_(0),
391    tiny_(1e-8),
392    veryTiny_(1e-20),
393    infty_(1e100),
394    exposeWarmStart_(false),
395    firstSolve_(true),
396    cutStrengthener_(NULL),
397    oaMessages_(),
398    oaHandler_(NULL),
399    newCutoffDecr(COIN_DBL_MAX)
400
401{
402   oaHandler_ = new OaMessageHandler;
403   oaHandler_->setLogLevel(0);
404}
405
406void 
407OsiTMINLPInterface::initialize(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
408                               Ipopt::SmartPtr<Ipopt::OptionsList> options,
409                               Ipopt::SmartPtr<Ipopt::Journalist> journalist,
410                               const std::string & prefix,
411                               Ipopt::SmartPtr<TMINLP> tminlp){
412  if(!IsValid(app_))
413     createApplication(roptions, options, journalist, prefix);
414  setModel(tminlp); 
415}
416
417void OsiTMINLPInterface::setSolver(Ipopt::SmartPtr<TNLPSolver> app){
418  app_ = app;
419  }
420
421
422void
423OsiTMINLPInterface::createApplication(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
424                                      Ipopt::SmartPtr<Ipopt::OptionsList> options,
425                                      Ipopt::SmartPtr<Ipopt::Journalist> journalist,
426                                      const std::string & prefix
427                                      )
428{
429  assert(!IsValid(app_));
430  int ival;
431  options->GetEnumValue("nlp_solver", ival, prefix.c_str());
432  Solver s = (Solver) ival;
433  if(s == EFilterSQP){
434    testOthers_ = false;;
435#ifdef COIN_HAS_FILTERSQP
436    app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
437#else
438   throw SimpleError("createApplication",
439                     "Bonmin not configured to run with FilterSQP.");
440#endif   
441
442   debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix)); 
443  }
444  else if(s == EIpopt){
445    testOthers_ = false;
446    app_ = new IpoptSolver(roptions, options, journalist, prefix);
447
448#ifdef COIN_HAS_FILTERSQP
449    debug_apps_.push_back(new Bonmin::FilterSolver(roptions, options, journalist, prefix));
450#endif
451  }
452  else if(s == EAll){
453#ifdef COIN_HAS_FILTERSQP
454    app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
455#else
456   throw SimpleError("createApplication",
457                     "Bonmin not configured to run with Ipopt.");
458#endif
459   debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix)); 
460    testOthers_ = true;
461  }
462  if (!app_->Initialize("")) {
463    throw CoinError("Error during initialization of app_","createApplication", "OsiTMINLPInterface");
464  }
465  for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin() ; 
466      i != debug_apps_.end() ; i++){
467    (*i)->Initialize("");
468  }
469  extractInterfaceParams();
470 
471}
472
473/** Facilitator to allocate a tminlp and an application. */
474void
475OsiTMINLPInterface::setModel(SmartPtr<TMINLP> tminlp)
476{
477  assert(IsValid(tminlp));
478  tminlp_ = tminlp;
479  problem_ = new TMINLP2TNLP(tminlp_);
480  feasibilityProblem_ = new TNLP2FPNLP
481        (SmartPtr<TNLP>(GetRawPtr(problem_)));
482  if(feasibility_mode_){
483    problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
484  }
485  else {
486    problem_to_optimize_ = GetRawPtr(problem_);
487  }
488}
489
490
491
492void
493OsiTMINLPInterface::readOptionFile(const std::string & fileName)
494{
495  if(IsValid(app_)){
496  std::ifstream is;
497  if (fileName != "") {
498    try {
499      is.open(fileName.c_str());
500    }
501    catch(std::bad_alloc) {
502      throw CoinError("Not enough memory to open option file.\n", "readOptionFile", "OsiTMINLPInterface");
503    }
504  }
505  options()->ReadFromStream(*app_->journalist(), is);
506  extractInterfaceParams();
507  }
508}
509
510/// Copy constructor
511OsiTMINLPInterface::OsiTMINLPInterface (const OsiTMINLPInterface &source):
512    OsiSolverInterface(source),
513    tminlp_(source.tminlp_),
514    problem_(NULL),
515    problem_to_optimize_(NULL),
516    feasibility_mode_(source.feasibility_mode_),
517    rowsense_(NULL),
518    rhs_(NULL),
519    rowrange_(NULL),
520    reducedCosts_(NULL),
521    OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
522    hasVarNamesFile_(source.hasVarNamesFile_),
523    nCallOptimizeTNLP_(0),
524    totalNlpSolveTime_(0),
525    totalIterations_(0),
526    maxRandomRadius_(source.maxRandomRadius_),
527    randomGenerationType_(source.randomGenerationType_),
528    max_perturbation_(source.max_perturbation_),
529    pushValue_(source.pushValue_),
530    numRetryInitial_(source.numRetryInitial_),
531    numRetryResolve_(source.numRetryResolve_),
532    numRetryInfeasibles_(source.numRetryInfeasibles_),
533    numRetryUnsolved_(source.numRetryUnsolved_),
534    dynamicCutOff_(source.dynamicCutOff_),
535    coeff_var_threshold_(source.coeff_var_threshold_),
536    first_perc_for_cutoff_decr_(source.first_perc_for_cutoff_decr_),
537    second_perc_for_cutoff_decr_(source.second_perc_for_cutoff_decr_),
538    messages_(),
539    pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
540    pretendSucceededNext_(source.pretendSucceededNext_),
541    hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
542    numIterationSuspect_(source.numIterationSuspect_),
543    hasBeenOptimized_(source.hasBeenOptimized_),
544    obj_(NULL),
545    feasibilityProblem_(NULL),
546    jRow_(NULL),
547    jCol_(NULL),
548    jValues_(NULL),
549    nnz_jac(source.nnz_jac),
550    constTypes_(NULL),
551//    constTypesNum_(NULL),
552    nNonLinear_(0),
553    tiny_(source.tiny_),
554    veryTiny_(source.veryTiny_),
555    infty_(source.infty_),
556    exposeWarmStart_(source.exposeWarmStart_),
557    firstSolve_(true),
558    cutStrengthener_(source.cutStrengthener_),
559    oaMessages_(),
560    oaHandler_(NULL),
561    newCutoffDecr(source.newCutoffDecr),
562    strong_branching_solver_(source.strong_branching_solver_)
563{
564   if(0 && defaultHandler()){
565     messageHandler()->setLogLevel(source.messageHandler()->logLevel());
566   }
567  //Pass in message handler
568  if(0 && source.messageHandler())
569    passInMessageHandler(source.messageHandler());
570   //Copy options from old application
571  if(IsValid(source.tminlp_)) {
572    problem_ = source.problem_->clone();
573    feasibilityProblem_ = new TNLP2FPNLP
574        (SmartPtr<TNLP>(GetRawPtr(problem_)), source.feasibilityProblem_);
575    if(feasibility_mode_)
576      problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
577    else
578      problem_to_optimize_ = GetRawPtr(problem_);
579    pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
580    pretendSucceededNext_ = source.pretendSucceededNext_;
581
582    setAuxiliaryInfo(source.getAuxiliaryInfo());
583    // Copy options from old application
584    app_ = source.app_->clone();
585    for(std::list<Ipopt::SmartPtr<TNLPSolver> >::const_iterator i = source.debug_apps_.begin();
586        i != source.debug_apps_.end() ; i++){
587      debug_apps_.push_back((*i)->clone());
588    }
589    testOthers_ = source.testOthers_;
590  }
591  else {
592    throw SimpleError("Don't know how to copy an empty IpoptInterface.",
593        "copy constructor");
594  }
595
596  warmstart_ = source.warmstart_ ? source.warmstart_->clone() : NULL;
597
598  if(source.obj_) {
599    obj_ = new double[source.getNumCols()];
600    CoinCopyN(source.obj_, source.getNumCols(), obj_);
601  }
602
603
604   oaHandler_ = new OaMessageHandler(*source.oaHandler_);;
605
606}
607
608OsiSolverInterface * 
609OsiTMINLPInterface::clone(bool copyData ) const
610{
611  if(copyData)
612    return new OsiTMINLPInterface(*this);
613  else return new OsiTMINLPInterface;
614}
615
616/// Assignment operator
617OsiTMINLPInterface & OsiTMINLPInterface::operator=(const OsiTMINLPInterface& rhs)
618{
619  if(this!= &rhs) {
620    OsiSolverInterface::operator=(rhs);
621    OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
622    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
623    totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
624    totalIterations_ = rhs.totalIterations_;
625    maxRandomRadius_ = rhs.maxRandomRadius_;
626    hasVarNamesFile_ = rhs.hasVarNamesFile_;
627    pushValue_ = rhs.pushValue_;
628
629    delete warmstart_;
630    warmstart_ = NULL;
631
632    if(IsValid(rhs.tminlp_)) {
633
634      tminlp_ = rhs.tminlp_;
635      problem_ = new TMINLP2TNLP(tminlp_);
636      problem_to_optimize_ = GetRawPtr(problem_);
637      app_ = rhs.app_->clone();
638
639      warmstart_ = rhs.warmstart_ ? rhs.warmstart_->clone() : NULL;
640
641      feasibilityProblem_ = new TNLP2FPNLP
642          (SmartPtr<TNLP>(GetRawPtr(problem_)));
643      nnz_jac = rhs.nnz_jac;
644
645      if(constTypes_ != NULL) {
646        delete [] constTypes_;
647        constTypes_ = NULL;
648      }
649      if(rhs.constTypes_ != NULL) {
650        constTypes_ = new TNLP::LinearityType[getNumRows()];
651        CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
652      }
653/*
654      if(constTypesNum_ != NULL) {
655        delete [] constTypesNum_;
656        constTypesNum_ = NULL;
657      }
658      if(rhs.constTypesNum_ != NULL) {
659        constTypesNum_ = new int[getNumRows()];
660        CoinCopyN(rhs.constTypesNum_, getNumRows(), constTypesNum_);
661      }
662*/
663      if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
664        jValues_ = new double [nnz_jac];
665        jCol_    = new Index [nnz_jac];
666        jRow_    = new Index [nnz_jac];
667        CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
668        CoinCopyN(rhs.jCol_    , nnz_jac,jCol_    );
669        CoinCopyN(rhs.jRow_    , nnz_jac,jRow_    );
670      }
671      else if(nnz_jac > 0) {
672        throw CoinError("Arrays for storing jacobian are inconsistant.",
673            "copy constructor",
674            "IpoptOAInterface");
675      }
676      tiny_ = rhs.tiny_;
677      veryTiny_ = rhs.veryTiny_;
678      infty_ = rhs.infty_;
679      exposeWarmStart_ = rhs.exposeWarmStart_;
680      newCutoffDecr = rhs.newCutoffDecr;
681
682    }
683    else {
684      tminlp_ =NULL;
685      problem_ = NULL;
686      app_ = NULL;
687      feasibilityProblem_ = NULL;
688    }
689
690
691    if(obj_) {
692      delete [] obj_;
693      obj_ = NULL;
694    }
695    if(rhs.obj_) {
696      obj_ = new double[rhs.getNumCols()];
697      CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
698    }
699
700    hasVarNamesFile_ = rhs.hasVarNamesFile_;
701
702    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
703    totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
704    totalIterations_ = rhs.totalIterations_;
705    maxRandomRadius_ = rhs.maxRandomRadius_;
706    pushValue_ = rhs.pushValue_;
707    numRetryInitial_ = rhs.numRetryInitial_;
708    numRetryResolve_ = rhs.numRetryResolve_;
709    numRetryInfeasibles_ = rhs.numRetryInfeasibles_;
710    numRetryUnsolved_ = rhs.numRetryUnsolved_;
711    pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
712    pretendSucceededNext_ = rhs.pretendSucceededNext_;
713    numIterationSuspect_ = rhs.numIterationSuspect_;
714
715    hasBeenOptimized_ = rhs.hasBeenOptimized_;
716    cutStrengthener_ = rhs.cutStrengthener_;
717
718    delete oaHandler_;
719    oaHandler_ = new OaMessageHandler(*rhs.oaHandler_);
720    strong_branching_solver_ = rhs.strong_branching_solver_;
721
722    dynamicCutOff_ = rhs.dynamicCutOff_;
723    coeff_var_threshold_ = rhs.coeff_var_threshold_;
724    first_perc_for_cutoff_decr_ = rhs.first_perc_for_cutoff_decr_;
725    second_perc_for_cutoff_decr_ = rhs.second_perc_for_cutoff_decr_;
726
727    freeCachedData();
728  }
729  return *this;
730}
731
732const SmartPtr<OptionsList> OsiTMINLPInterface::options() const
733{
734  if(!IsValid(app_)) {
735    messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
736    return NULL;
737  }
738  else
739    return app_->options();
740}
741
742SmartPtr<OptionsList> OsiTMINLPInterface::options()
743{
744  if(!IsValid(app_)) {
745    messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
746    return NULL;
747  }
748  else
749    return app_->options();
750}
751
752/// Destructor
753OsiTMINLPInterface::~OsiTMINLPInterface ()
754{
755  freeCachedData();
756  delete [] jRow_;
757  delete [] jCol_;
758  delete [] jValues_;
759  delete [] constTypes_;
760  delete [] obj_;
761  delete oaHandler_;
762  delete warmstart_;
763}
764
765void
766OsiTMINLPInterface::freeCachedColRim()
767{
768    if(reducedCosts_!=NULL) {
769    delete []  reducedCosts_;
770    reducedCosts_ = NULL;
771  }
772
773}
774
775void
776OsiTMINLPInterface::freeCachedRowRim()
777{
778  if(rowsense_!=NULL) {
779    delete []  rowsense_;
780    rowsense_ = NULL;
781  }
782  if(rhs_!=NULL) {
783    delete []  rhs_;
784    rhs_ = NULL;
785  }
786  if(rowrange_!=NULL) {
787    delete []  rowrange_;
788    rowrange_ = NULL;
789  }
790  //   if(dualsol_!=NULL)
791  //     {
792  //       delete []  dualsol_; dualsol_ = NULL;
793  //     }
794}
795
796void
797OsiTMINLPInterface::freeCachedData()
798{
799  freeCachedColRim();
800  freeCachedRowRim();
801}
802
803const char * OsiTMINLPInterface::OPT_SYMB="OPT";
804const char * OsiTMINLPInterface::FAILED_SYMB="FAILED";
805const char * OsiTMINLPInterface::UNBOUND_SYMB="UNBOUNDED";
806const char * OsiTMINLPInterface::INFEAS_SYMB="INFEAS";
807const char * OsiTMINLPInterface::TIME_SYMB="TIME";
808///////////////////////////////////////////////////////////////////
809// WarmStart Information                                                                           //
810///////////////////////////////////////////////////////////////////
811
812
813void
814OsiTMINLPInterface::resolveForCost(int numsolve, bool keepWarmStart)
815{
816  // This method assumes that a problem has just been solved and we try for a
817  // different solution. So disregard (in fact, clear out) any warmstart info
818  // we might have, and acquire a new one before returning.
819  delete warmstart_;
820  warmstart_ = NULL;
821 
822  double * of_current = (numsolve > 0) ? new double[numsolve]: NULL;
823  int num_failed, num_infeas;
824  double mean, std_dev, var_coeff;
825  double min = DBL_MAX;
826  double max = -DBL_MAX;
827
828  Coin::SmartPtr<SimpleReferencedPtr<CoinWarmStart> > ws_backup = NULL;
829  if(!exposeWarmStart_ && keepWarmStart){
830    //if warm start is not exposed, we need to store the current starting point to
831    //restore it at the end of the method
832    ws_backup = make_referenced(app_->getUsedWarmStart(problem_));
833  }
834  /** Save current optimum. */
835  vector<double> point(getNumCols()*3+ getNumRows());
836  double bestBound = (isProvenOptimal()) ? getObjValue() : DBL_MAX;
837  CoinCopyN(getColSolution(),
838      getNumCols(), point());
839  CoinCopyN(getRowPrice(),
840      2*getNumCols()+ getNumRows(),
841      point() + getNumCols());
842  TNLPSolver::ReturnStatus savedStatus = optimizationStatus_;
843  if(isProvenOptimal())
844    messageHandler()->message(SOLUTION_FOUND,
845        messages_)
846    <<1<<getObjValue()<<bestBound
847    <<CoinMessageEol;
848  else
849    messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
850        messages_)
851    <<1
852    <<CoinMessageEol;
853
854  num_failed = 0;
855  num_infeas = 0;
856  mean = 0;
857
858  for(int f = 0; f < numsolve ; f++) {
859    messageHandler()->message(WARNING_RESOLVING,
860        messages_)
861    <<f+1<< CoinMessageEol ;
862    randomStartingPoint();
863    solveAndCheckErrors(0,0,"resolve cost");
864
865
866    char c=' ';
867    //Is solution better than previous
868    if(isProvenOptimal() &&
869        getObjValue()<bestBound) {
870      c='*';
871      messageHandler()->message(BETTER_SOL, messages_)<<getObjValue()<<f+1<< CoinMessageEol;
872      CoinCopyN(getColSolution(),
873          getNumCols(), point());
874      CoinCopyN(getRowPrice(),
875          2*getNumCols()+ getNumRows(),
876          point() + getNumCols());
877      bestBound = getObjValue();
878      savedStatus = optimizationStatus_;
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
946  if(dynamicCutOff_)
947  {
948     if(var_coeff<0.1)
949     {
950        setNewCutoffDecr(mean*first_perc_for_cutoff_decr_);
951     }
952     else
953     {
954        setNewCutoffDecr(mean*second_perc_for_cutoff_decr_);
955     }
956  }
957     
958
959  problem_->Set_x_sol(getNumCols(),point());
960  problem_->Set_dual_sol(point.size()-getNumCols(), point() + getNumCols());
961  problem_->set_obj_value(bestBound);
962  optimizationStatus_ = savedStatus;
963  hasBeenOptimized_ = true;
964
965  if(!exposeWarmStart_ && keepWarmStart) {
966    app_->setWarmStart(ws_backup->ptr(), problem_);
967  }
968}
969
970void
971OsiTMINLPInterface::resolveForRobustness(int numsolve)
972{
973  // This method assumes that a problem has just been solved and we try for a
974  // different solution. So disregard (in fact, clear out) any warmstart info
975  // we might have, and acquire a new one before returning.
976  delete warmstart_;
977  warmstart_ = NULL;
978 
979
980  CoinWarmStart * ws_backup = NULL;
981  if(!exposeWarmStart_){
982    //if warm start is not exposed, we need to store the current starting point to
983    //restore it at the end of the method
984    ws_backup = app_->getUsedWarmStart(problem_);
985  }
986  //std::cerr<<"Resolving the problem for robustness"<<std::endl;
987  //First remove warm start point and resolve
988  app_->disableWarmStart();
989  problem()->resetStartingPoint();
990  messageHandler()->message(WARNING_RESOLVING,
991      messages_)
992  <<1<< CoinMessageEol ;
993  solveAndCheckErrors(0,0,"resolve robustness");
994
995
996  char c='*';
997  if(isAbandoned()) {
998    c=' ';
999  }
1000  messageHandler()->message(LOG_LINE, messages_)
1001  <<c<<1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<
1002    app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
1003
1004
1005  if(!isAbandoned()) {
1006    messageHandler()->message(WARN_SUCCESS_WS, messages_) << CoinMessageEol ;
1007    // re-enable warmstart and get it
1008    app_->enableWarmStart();
1009    if (!exposeWarmStart_) {
1010      app_->setWarmStart(ws_backup, problem_);
1011      delete ws_backup;
1012    }
1013    return; //we won go on
1014  }
1015
1016  //still unsolved try again with different random starting points
1017  for(int f = 0; f < numsolve ; f++) {
1018    messageHandler()->message(WARNING_RESOLVING,
1019        messages_)
1020    <<f+2<< CoinMessageEol ;
1021
1022    randomStartingPoint();
1023    solveAndCheckErrors(0,0,"resolve robustness");
1024
1025
1026    messageHandler()->message(IPOPT_SUMMARY, messages_)
1027    <<"resolveForRobustness"<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
1028
1029
1030    char c='*';
1031    if(isAbandoned()) {
1032      c=' ';
1033    }
1034    messageHandler()->message(LOG_LINE, messages_)
1035    <<c<<f+2<<statusAsString()<<getObjValue()
1036    <<app_->IterationCount()<<app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
1037
1038
1039    if(!isAbandoned()) {
1040      messageHandler()->message(WARN_SUCCESS_RANDOM, messages_)
1041        << f+2 << CoinMessageEol ;
1042      // re-enable warmstart and get it
1043      app_->enableWarmStart();
1044      if (!exposeWarmStart_) {
1045        app_->setWarmStart(ws_backup, problem_);
1046        delete ws_backup;
1047      }
1048       
1049      return; //we have found a solution and continue
1050    }
1051  }
1052
1053
1054  if(!exposeWarmStart_){
1055    app_->setWarmStart(ws_backup, problem_);
1056    delete ws_backup;
1057  }
1058  if(pretendFailIsInfeasible_) {
1059    pretendSucceededNext_ = true;
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  for(int i = 0 ; i < m ; i++) {
2354    if(constTypes_[i] == TNLP::NON_LINEAR) {
2355      //If constraint is range not binding prepare to remove it
2356      if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty && fabs(duals[i]) == 0.)
2357      {
2358        nonBindings[numNonBindings++] = i;
2359        continue;
2360      }
2361      else
2362        if(rowLower[i] > - nlp_infty){
2363          rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
2364          if(! WarnedForNonConvexOa && rowUpper[i] < nlp_infty){
2365             messageHandler()->message(WARNING_NON_CONVEX_OA, messages_)<<CoinMessageEol;
2366             WarnedForNonConvexOa = true;
2367          }
2368        }
2369      else
2370        rowLow[i] = - infty;
2371      if(rowUpper[i] < nlp_infty)
2372        rowUp[i] =  (rowUpper[i] - g[i]) + 1e-07;
2373      else
2374        rowUp[i] = infty;
2375     
2376      //If equality or ranged constraint only add one side by looking at sign of dual multiplier
2377      if(rowLower[i] > -nlp_infty && rowUpper[i] < nlp_infty)
2378      {
2379        if(duals[i] >= 0.)// <= inequality
2380          rowLow[i] = - infty;
2381        if(duals[i] <= 0.)// >= inequality
2382          rowUp[i] = infty;
2383      }
2384    }
2385    else {
2386      if(rowLower[i] > -nlp_infty){
2387      //   printf("Lower %g ", rowLower[i]);
2388         //rowLow[i] = (rowLower[i] - g[i]);
2389         rowLow[i] = (rowLower[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]);
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          if(rowLow[jRow_[i]] > - infty)
2415          rowLow[jRow_[i]] += jValues_[i] * x[jCol_ [i]];
2416          if(rowUp[jRow_[i]] < infty)
2417          rowUp[jRow_[i]] += jValues_[i] *x[jCol_[i]];
2418       }
2419    }
2420    else {
2421      //double value = jValues_[i] * x[jCol_[i]];
2422      //rowLow[jRow_[i]] += value;
2423      //rowUp[jRow_[i]] += value;
2424    } 
2425  }
2426  CoinPackedMatrix mat(true, jRow_, jCol_, jValues_, nnz_jac_g);
2427  mat.setDimensions(m,n); // In case matrix was empty, this should be enough
2428 
2429  //remove non-bindings equality constraints
2430  mat.deleteRows(numNonBindings, nonBindings());
2431
2432  int numcols=getNumCols();
2433  vector<double> obj(numcols);
2434  for(int i = 0 ; i < numcols ; i++)
2435    obj[i] = 0.;
2436 
2437 
2438  si.loadProblem(mat, getColLower(), getColUpper(), obj(), rowLow(), rowUp());
2439  for(int i = 0 ; i < getNumCols() ; i++) {
2440    if(isInteger(i))
2441      si.setInteger(i);
2442  }
2443  if(getObj) {
2444     bool addObjVar = false;
2445     if(problem_->hasLinearObjective()){
2446       double zero;
2447       vector<double> x0(n,0.);
2448       problem_to_optimize_->eval_f(n, x0(), 1, zero);
2449       si.setDblParam(OsiObjOffset, -zero);
2450       //Copy the linear objective and don't create a dummy variable.
2451       problem_to_optimize_->eval_grad_f(n, x, 1,obj());
2452       si.setObjective(obj());
2453    }
2454    else {
2455      addObjVar = true;
2456   }
2457   
2458   if(addObjVar){
2459      addObjectiveFunction(si, x);
2460    }
2461  }
2462
2463}
2464
2465void 
2466OsiTMINLPInterface::addObjectiveFunction(OsiSolverInterface &si, 
2467                                         const double *x){
2468  const double * colLower = getColLower();
2469  const double * colUpper = getColUpper();
2470  int numcols = getNumCols();
2471  assert(numcols == si.getNumCols() );
2472  vector<double> obj(numcols);
2473      problem_to_optimize_->eval_grad_f(numcols, x, 1,obj());
2474      //add variable alpha
2475      //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
2476      CoinPackedVector a;
2477      si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
2478 
2479      // Now get the objective cuts
2480      // get the gradient, pack it and add the cut
2481      double ub;
2482      problem_to_optimize_->eval_f(numcols, x, 1, ub);
2483      ub*=-1;
2484      double lb = -1e300;
2485      CoinPackedVector objCut;
2486      CoinPackedVector * v = &objCut;
2487      v->reserve(numcols+1);
2488      for(int i = 0; i<numcols ; i++) {
2489       if(si.getNumRows())
2490       {
2491        if(cleanNnz(obj[i],colLower[i], colUpper[i],
2492            -getInfinity(), 0,
2493            x[i],
2494            lb,
2495            ub, tiny_, veryTiny_)) {
2496          v->insert(i,obj[i]);
2497          lb += obj[i] * x[i];
2498          ub += obj[i] * x[i];
2499        }
2500       }
2501       else //Unconstrained problem can not put clean coefficient
2502       {
2503           if(cleanNnz(obj[i],colLower[i], colUpper[i],
2504            -getInfinity(), 0,
2505            x[i],
2506            lb,
2507            ub, 1e-03, 1e-08)) {
2508          v->insert(i,obj[i]);
2509          lb += obj[i] * x[i];
2510          ub += obj[i] * x[i];
2511           }
2512       }
2513      }
2514    v->insert(numcols,-1);
2515    si.addRow(objCut, lb, ub);
2516    }
2517
2518/** Add a collection of linear cuts to problem formulation.*/
2519void 
2520OsiTMINLPInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
2521{ 
2522  if(numberCuts)
2523    freeCachedRowRim();
2524  const OsiRowCut ** cutsPtrs = new const OsiRowCut*[numberCuts];
2525  for(int i = 0 ; i < numberCuts ; i++)
2526  {
2527    cutsPtrs[i] = &cuts[i];
2528  }
2529  problem_->addCuts(numberCuts, cutsPtrs);
2530  delete [] cutsPtrs;
2531}
2532
2533void
2534OsiTMINLPInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
2535    const char * whereFrom)
2536{
2537  totalNlpSolveTime_-=CoinCpuTime();
2538  if(warmStarted)
2539    optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
2540  else
2541    optimizationStatus_ = app_->OptimizeTNLP(GetRawPtr(problem_to_optimize_));
2542  totalNlpSolveTime_+=CoinCpuTime();
2543  nCallOptimizeTNLP_++;
2544  hasBeenOptimized_ = true;
2545 
2546 
2547  //Options should have been printed if not done already turn off Ipopt output
2548  if(!hasPrintedOptions) {
2549    hasPrintedOptions = 1;
2550    //app_->Options()->SetIntegerValue("print_level",0, true, true);
2551    app_->options()->SetStringValue("print_user_options","no", false, true);
2552  }
2553 
2554  bool otherDisagree = false ;
2555#if 0
2556  if(optimizationStatus_ == TNLPSolver::notEnoughFreedom)//Too few degrees of freedom
2557  {
2558    (*messageHandler())<<"Too few degrees of freedom...."<<CoinMessageEol;
2559    int numrows = getNumRows();
2560    int numcols = getNumCols();
2561   
2562    const double * colLower = getColLower();
2563    const double * colUpper = getColUpper();
2564   
2565   
2566    const double * rowLower = getRowLower();
2567    const double * rowUpper = getRowUpper();
2568   
2569    int numberFixed = 0;
2570    for(int i = 0 ; i < numcols ; i++)
2571    {
2572      if(colUpper[i] - colLower[i] <= INT_BIAS)
2573            {
2574              numberFixed++;
2575            }
2576    }
2577    int numberEqualities = 0;
2578    for(int i = 0 ; i < numrows ; i++)
2579    {
2580      if(rowUpper[i] - rowLower[i] <= INT_BIAS)
2581            {
2582              numberEqualities++;
2583            }     
2584    }
2585    if(numcols - numberFixed > numberEqualities || numcols < numberEqualities)
2586    {
2587      std::string probName;
2588      getStrParam(OsiProbName, probName);
2589      throw newUnsolvedError(app_->errorCode(), problem_, probName);
2590    }
2591    double * saveColLow = CoinCopyOfArray(getColLower(), getNumCols());
2592    double * saveColUp = CoinCopyOfArray(getColUpper(), getNumCols());
2593   
2594    for(int i = 0; i < numcols && numcols - numberFixed <= numberEqualities ; i++)
2595    {
2596      if(colUpper[i] - colLower[i] <= INT_BIAS)
2597            {
2598              setColLower(i, saveColLow[i]-1e-06);
2599              setColUpper(i, saveColLow[i]+1e-06);
2600              numberFixed--;
2601            }
2602    }
2603    solveAndCheckErrors(warmStarted, throwOnFailure, whereFrom);
2604    //restore
2605    for(int i = 0; i < numcols && numcols - numberFixed < numrows ; i++)
2606    {
2607      problem_->SetVariableLowerBound(i,saveColLow[i]);
2608      problem_->SetVariableUpperBound(i,saveColUp[i]);
2609    }
2610    delete [] saveColLow;
2611    delete [] saveColUp;
2612    return;
2613  }
2614  else
2615#endif
2616    if(!app_->isRecoverable(optimizationStatus_))//Solver failed and the error can not be recovered, throw it
2617    {
2618      std::string probName;
2619      getStrParam(OsiProbName, probName);
2620      throw newUnsolvedError(app_->errorCode(), problem_, probName);
2621    }
2622    else if(testOthers_ && !app_->isError(optimizationStatus_)){
2623      Ipopt::SmartPtr<TMINLP2TNLP> problem_copy = problem_->clone();
2624      //Try other solvers and see if they agree
2625      int f =1;
2626      for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2627          i != debug_apps_.end() ; i++){
2628        TNLPSolver::ReturnStatus otherStatus = (*i)->OptimizeTNLP(GetRawPtr(problem_copy));
2629       messageHandler()->message(LOG_LINE, messages_)
2630         <<'d'<<f++<<statusAsString(otherStatus)<<problem_copy->obj_value()
2631         <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
2632        if(!(*i)->isError(otherStatus)){
2633           CoinRelFltEq eq(1e-05);
2634           if(otherStatus != optimizationStatus_){
2635             otherDisagree = true;
2636             messageHandler()->message(SOLVER_DISAGREE_STATUS, messages_)
2637             <<app_->solverName()<<statusAsString()
2638             <<(*i)->solverName()<<statusAsString(otherStatus)<<CoinMessageEol; 
2639           }
2640           else if(isProvenOptimal() && !eq(problem_->obj_value(),problem_copy->obj_value()))
2641           {
2642             otherDisagree = true;
2643             messageHandler()->message(SOLVER_DISAGREE_VALUE, messages_)
2644             <<app_->solverName()<<problem_->obj_value()
2645             <<(*i)->solverName()<<problem_copy->obj_value()<<CoinMessageEol; 
2646           }
2647        }
2648     }
2649  }
2650  else if(app_->isError(optimizationStatus_) && ! debug_apps_.empty()){
2651      int f =1;
2652      for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2653          i != debug_apps_.end() && app_->isError(optimizationStatus_) ; i++){
2654        optimizationStatus_ = (*i)->OptimizeTNLP(GetRawPtr(problem_));
2655        messageHandler()->message(LOG_LINE, messages_)
2656          <<'d'<<f++<<statusAsString(optimizationStatus_)<<problem_->obj_value()
2657          <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
2658      }
2659  }
2660  try{
2661    totalIterations_ += app_->IterationCount();
2662  }
2663  catch(SimpleError &E)
2664  {
2665    if (throwOnFailure)//something failed throw
2666    {
2667      throw SimpleError("No statistics available from Ipopt",whereFrom);
2668    }
2669    else {
2670      return;
2671    }
2672  }
2673  if(problem_->hasUpperBoundingObjective()){//Check if solution is integer and recompute objective value using alternative objective function
2674    const double * sol = getColSolution();
2675    bool integerSol = true;
2676    double intTol = 1e-08;
2677    if(objects()){
2678      int nObjects = numberObjects();
2679      OsiObject ** object = objects();
2680      for(int i = 0 ; i< nObjects ; i++){
2681        int dummy;
2682        if(object[i]->infeasibility(this,dummy) > intTol)
2683        {
2684          integerSol=false;
2685          break;
2686        }
2687      }
2688    }
2689    else{//Only works for integer constraints
2690      int numcols = getNumCols();
2691      for(int i = 0 ; i < numcols ; i++){
2692        if(isInteger(i) || isBinary(i)){
2693          if(fabs(sol[i] - floor(sol[i]+0.5)) > intTol){
2694            integerSol = false;
2695            break;
2696          }
2697        }
2698      }
2699    }
2700    if(integerSol&&isProvenOptimal()){
2701      double help= problem_->evaluateUpperBoundingFunction(sol);
2702     
2703
2704      OsiAuxInfo * auxInfo = getAuxiliaryInfo();
2705      Bonmin::AuxInfo * bonInfo = dynamic_cast<Bonmin::AuxInfo *>(auxInfo);
2706      if(bonInfo!=0)
2707      {
2708       
2709        if(help<bonInfo->bestObj2())
2710        {
2711          bonInfo->setBestObj2(help);
2712          bonInfo->setBestSolution2(getNumCols(), const_cast<double *>(getColSolution()));
2713
2714           messageHandler()->message(ALTERNATE_OBJECTIVE, messages_)
2715           <<help<<CoinMessageEol;
2716        }
2717      }
2718      else {
2719        printf("\nWARNING: the algorithm selected does not consider the second objective function\n");
2720      }
2721    }
2722  }
2723  messageHandler()->message(IPOPT_SUMMARY, messages_)
2724    <<whereFrom<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
2725 
2726  if((nCallOptimizeTNLP_ % 20) == 1)
2727    messageHandler()->message(LOG_HEAD, messages_)<<CoinMessageEol;
2728 
2729 
2730  if ( (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) ||
2731       ( otherDisagree )){
2732    messageHandler()->message(SUSPECT_PROBLEM,
2733                              messages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
2734    std::string subProbName;
2735    getStrParam(OsiProbName, subProbName);
2736    std::ostringstream os;
2737    os<<"_"<<nCallOptimizeTNLP_;
2738    subProbName+=os.str();
2739    problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
2740  }
2741 
2742}
2743
2744////////////////////////////////////////////////////////////////////
2745// Solve Methods                                                  //
2746////////////////////////////////////////////////////////////////////
2747void OsiTMINLPInterface::initialSolve()
2748{
2749   initialSolve("");
2750}
2751
2752////////////////////////////////////////////////////////////////////
2753// Solve Methods                                                  //
2754////////////////////////////////////////////////////////////////////
2755void OsiTMINLPInterface::initialSolve(const char * whereFrom)
2756{
2757  assert(IsValid(app_));
2758  assert(IsValid(problem_));
2759
2760  // Discard warmstart_ if we had one
2761  delete warmstart_;
2762  warmstart_ = NULL;
2763 
2764  if(!hasPrintedOptions) {
2765    int printOptions;
2766    app_->options()->GetEnumValue("print_user_options",printOptions,app_->prefix());
2767    if(printOptions)
2768      app_->options()->SetStringValue("print_user_options","yes",true,true);
2769  }
2770  if(exposeWarmStart_)
2771    app_->disableWarmStart(); 
2772  solveAndCheckErrors(0,1,"initialSolve");
2773 
2774  //Options should have been printed if not done already turn off Ipopt output
2775  if(!hasPrintedOptions) {
2776    hasPrintedOptions = 1;
2777    app_->options()->SetStringValue("print_user_options","no");
2778    app_->options()->SetIntegerValue("print_level",0);
2779  }
2780 
2781  messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
2782                                                      <<statusAsString()
2783                                                      <<getObjValue()
2784                                                      <<app_->IterationCount()
2785                                                      <<app_->CPUTime()
2786                                                      <<whereFrom<<CoinMessageEol;
2787 
2788  int numRetry = firstSolve_ ? numRetryInitial_ : numRetryResolve_;
2789  if(isAbandoned()) {
2790    resolveForRobustness(numRetryUnsolved_);
2791  }
2792  else if(numRetry)
2793    {
2794      resolveForCost(numRetry, numRetryInitial_ > 0);
2795      /** Only do it once at the root.*/
2796      numRetryInitial_ = 0;
2797    }
2798  firstSolve_ = false;
2799
2800  // if warmstart_ is not empty then had to use resolveFor... and that created
2801  // the warmstart at the end, and we have nothing to do here. Otherwise...
2802  if (! warmstart_ && ! isAbandoned()) {
2803    if (exposeWarmStart_) {
2804      warmstart_ = app_->getWarmStart(problem_);
2805    }
2806  }
2807}
2808
2809/** Resolve the continuous relaxation after problem modification. */
2810void
2811OsiTMINLPInterface::resolve(){
2812   resolve("");
2813}
2814/** Resolve the continuous relaxation after problem modification.
2815 * \note for Ipopt, same as resolve */
2816void
2817OsiTMINLPInterface::resolve(const char * whereFrom)
2818{
2819  assert(IsValid(app_));
2820  assert(IsValid(problem_));
2821  int has_warmstart = warmstart_ == NULL ? 0 : 1;
2822  if(warmstart_ == NULL) has_warmstart = 0;
2823  else if(!app_->warmStartIsValid(warmstart_)) has_warmstart = 1;
2824  else has_warmstart = 2;
2825  if (has_warmstart < 2) {
2826    // empty (or unrecognized) warmstart
2827    initialSolve(whereFrom);
2828    return;
2829  }
2830  app_->setWarmStart(warmstart_, problem_);
2831  delete warmstart_;
2832  warmstart_ = NULL;
2833
2834  if (INT_BIAS > 0.) {
2835    app_->options()->SetStringValue("warm_start_same_structure", "yes");
2836  }
2837  else {
2838    app_->options()->SetStringValue("warm_start_same_structure", "no");
2839  }
2840
2841  if(problem_->duals_init() != NULL)
2842    app_->enableWarmStart();
2843  else app_->disableWarmStart();
2844  solveAndCheckErrors(1,1,"resolve");
2845 
2846  messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
2847                                                <<statusAsString()
2848                                                <<getObjValue()
2849                                                <<app_->IterationCount()
2850                                                <<app_->CPUTime()
2851                                                <<whereFrom
2852                                                <<"totot"
2853                                                <<CoinMessageEol;
2854 
2855  if(isAbandoned()) {
2856    resolveForRobustness(numRetryUnsolved_);
2857  }
2858  else if(numRetryResolve_ ||
2859          (numRetryInfeasibles_ && isProvenPrimalInfeasible() ))
2860    resolveForCost(std::max(numRetryResolve_, numRetryInfeasibles_), 0);
2861
2862  // if warmstart_ is not empty then had to use resolveFor... and that created
2863  // the warmstart at the end, and we have nothing to do here. Otherwise...
2864  if (! warmstart_ && ! isAbandoned()) {
2865    if (exposeWarmStart_) {
2866      warmstart_ = app_->getWarmStart(problem_);
2867    }
2868  }
2869}
2870
2871
2872////////////////////////////////////////////////////////////////
2873// Methods returning info on how the solution process terminated  //
2874////////////////////////////////////////////////////////////////
2875/// Are there a numerical difficulties?
2876bool OsiTMINLPInterface::isAbandoned() const
2877{
2878  if (pretendSucceededNext_)
2879    return false;
2880  pretendSucceededNext_ = false;
2881  return ( 
2882        (optimizationStatus_==TNLPSolver::iterationLimit)||
2883        (optimizationStatus_==TNLPSolver::computationError)||
2884        (optimizationStatus_==TNLPSolver::illDefinedProblem)||
2885        (optimizationStatus_==TNLPSolver::illegalOption)||
2886        (optimizationStatus_==TNLPSolver::externalException)|
2887        (optimizationStatus_==TNLPSolver::exception)
2888      );
2889}
2890
2891/// Is optimality proven?
2892bool OsiTMINLPInterface::isProvenOptimal() const
2893{
2894  return (optimizationStatus_==TNLPSolver::solvedOptimal) ||
2895          (optimizationStatus_==TNLPSolver::solvedOptimalTol);
2896}
2897/// Is primal infeasiblity proven?
2898bool OsiTMINLPInterface::isProvenPrimalInfeasible() const
2899{
2900  return (optimizationStatus_ == TNLPSolver::provenInfeasible);
2901}
2902/// Is dual infeasiblity proven?
2903bool OsiTMINLPInterface::isProvenDualInfeasible() const
2904{
2905  return (optimizationStatus_ == TNLPSolver::unbounded);
2906}
2907/// Is the given primal objective limit reached?
2908bool OsiTMINLPInterface::isPrimalObjectiveLimitReached() const
2909{
2910  (*handler_)<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2911  return 0;
2912}
2913/// Is the given dual objective limit reached?
2914bool OsiTMINLPInterface::isDualObjectiveLimitReached() const
2915{
2916  //  (*messageHandler_)<<"Warning : isDualObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2917  return (optimizationStatus_==TNLPSolver::unbounded);
2918
2919}
2920/// Iteration limit reached?
2921bool OsiTMINLPInterface::isIterationLimitReached() const
2922{
2923  return (optimizationStatus_==TNLPSolver::iterationLimit);
2924}
2925
2926void
2927OsiTMINLPInterface::extractInterfaceParams()
2928{
2929  if (IsValid(app_)) {
2930    int logLevel;
2931    app_->options()->GetIntegerValue("nlp_log_level", logLevel,app_->prefix());
2932    messageHandler()->setLogLevel(logLevel);
2933
2934#ifdef COIN_HAS_FILTERSQP
2935    FilterSolver * filter = dynamic_cast<FilterSolver *>(GetRawPtr(app_));
2936
2937    bool is_given =
2938#endif
2939      app_->options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,app_->prefix());
2940
2941#ifdef COIN_HAS_FILTERSQP
2942    if(filter && !is_given){
2943      // overwriting default for filter
2944      maxRandomRadius_ = 10.;
2945    }
2946#endif
2947   
2948   int oaCgLogLevel = 0;
2949   app_->options()->GetIntegerValue("oa_cuts_log_level", oaCgLogLevel,app_->prefix());
2950   oaHandler_->setLogLevel(oaCgLogLevel); 
2951       
2952    int exposeWs = false;
2953    app_->options()->GetEnumValue("warm_start", exposeWs, app_->prefix());
2954    setExposeWarmStart(exposeWs > 0);
2955   
2956    app_->options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,app_->prefix());
2957    app_->options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,app_->prefix());
2958    app_->options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,app_->prefix());
2959    app_->options()->GetIntegerValue("num_resolve_at_infeasibles", numRetryInfeasibles_,app_->prefix());
2960    app_->options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,app_->prefix());
2961    app_->options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,app_->prefix());
2962    app_->options()->GetNumericValue
2963    ("warm_start_bound_frac" ,pushValue_,app_->prefix());
2964    app_->options()->GetNumericValue("tiny_element",tiny_,app_->prefix());
2965    app_->options()->GetNumericValue("very_tiny_element",veryTiny_,app_->prefix());
2966    app_->options()->GetNumericValue("random_point_perturbation_interval",max_perturbation_,app_->prefix());
2967    app_->options()->GetEnumValue("random_point_type",randomGenerationType_,app_->prefix());
2968    int cut_strengthening_type;
2969    app_->options()->GetEnumValue("cut_strengthening_type", cut_strengthening_type,app_->prefix());
2970
2971    if (cut_strengthening_type != CS_None) {
2972      // TNLP solver to be used in the cut strengthener
2973      cutStrengthener_ = new CutStrengthener(app_->clone(), app_->options());
2974    }
2975  }
2976}
2977
2978void
2979OsiTMINLPInterface::SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)
2980{
2981  strong_branching_solver_ = strong_branching_solver;
2982}
2983
2984  //#define STRONG_COMPARE
2985#ifdef STRONG_COMPARE
2986  static double objorig;
2987#endif
2988
2989void
2990OsiTMINLPInterface::markHotStart()
2991{
2992  if (IsValid(strong_branching_solver_)) {
2993#ifdef STRONG_COMPARE
2994    // AWDEBUG
2995    OsiSolverInterface::markHotStart();
2996    objorig = getObjValue();
2997#endif
2998    optimizationStatusBeforeHotStart_ = optimizationStatus_;
2999    strong_branching_solver_->markHotStart(this);
3000  }
3001  else {
3002    // Default Implementation
3003    OsiSolverInterface::markHotStart();
3004  }
3005}
3006
3007void
3008OsiTMINLPInterface::solveFromHotStart()
3009{
3010  if (IsValid(strong_branching_solver_)) {
3011#ifdef STRONG_COMPARE
3012    // AWDEBUG
3013    OsiSolverInterface::solveFromHotStart();
3014    double obj_nlp = getObjValue() - objorig;
3015#endif
3016    optimizationStatus_ = strong_branching_solver_->solveFromHotStart(this);
3017    hasBeenOptimized_ = true;
3018#ifdef STRONG_COMPARE
3019    double obj_other = getObjValue() - objorig;
3020    printf("AWDEBUG: Strong Branching results: NLP = %15.8e Other = %15.8e\n",
3021           obj_nlp, obj_other);
3022#endif
3023  }
3024  else {
3025    // Default Implementation
3026    OsiSolverInterface::solveFromHotStart();
3027  }
3028}
3029
3030void
3031OsiTMINLPInterface::unmarkHotStart()
3032{
3033  if (IsValid(strong_branching_solver_)) {
3034#ifdef STRONG_COMPARE
3035    // AWDEBUG
3036    OsiSolverInterface::unmarkHotStart();
3037#endif
3038    strong_branching_solver_->unmarkHotStart(this);
3039    optimizationStatus_ = optimizationStatusBeforeHotStart_;
3040  }
3041  else {
3042    // Default Implementation
3043    OsiSolverInterface::unmarkHotStart();
3044  }
3045}
3046
3047const double * OsiTMINLPInterface::getObjCoefficients() const
3048{
3049  const int n = getNumCols();
3050  delete [] obj_;
3051  obj_ = NULL;
3052  obj_ = new double[n];
3053
3054  bool new_x = true;
3055  const double* x_sol = problem_->x_sol();
3056  bool retval = problem_->eval_grad_f(n, x_sol, new_x, obj_);
3057 
3058  if (!retval) {
3059    // Let's see if that happens - it will cause a crash
3060    fprintf(stderr, "ERROR WHILE EVALUATING GRAD_F in OsiTMINLPInterface::getObjCoefficients()\n");
3061    delete [] obj_;
3062    obj_ = NULL;
3063  }
3064
3065  return obj_;
3066}
3067
3068  /** Sets the TMINLP2TNLP to be used by the interface.*/
3069  void 
3070  OsiTMINLPInterface::use(Ipopt::SmartPtr<TMINLP2TNLP> tminlp2tnlp){
3071     problem_ = tminlp2tnlp;
3072     problem_to_optimize_ = GetRawPtr(problem_);
3073     feasibilityProblem_->use(GetRawPtr(tminlp2tnlp));}
3074
3075}/** end namespace Bonmin*/
3076
Note: See TracBrowser for help on using the repository browser.