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

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

Remove some unused option in OA add Claudi d'Ambrosio code

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