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

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

Small documenting fix

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