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

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

merge chgset 1730 from stable/1.4: some fixes to options documentation

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