source: branches/devel/Bonmin/src/IpoptInterface/IpoptInterface.cpp @ 26

Last change on this file since 26 was 26, checked in by andreasw, 13 years ago

Fix ticket #2, improve behavior of OA in non-convex problems

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 73.6 KB
Line 
1// (C) Copyright International Business Machines Corporation and Carnegie Mellon University 2004
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// Pierre Bonami, Carnegie Mellon University,
7// Carl D. Laird, Carnegie Mellon University,
8// Andreas Waechter, International Business Machines Corporation
9//
10// Date : 12/01/2004
11
12
13#include "IpoptInterface.hpp"
14#include "TMINLP.hpp"
15#include "IpCbcColReader.hpp"
16#include "CoinTime.hpp"
17#include "IpoptWarmStart.hpp"
18#include "IpoptIntegerBias.hpp"
19#include "IpoptInteriorWarmStarter.hpp"
20#include <string>
21#include <sstream>
22
23#include "IpSolveStatistics.hpp"
24
25//For the branch and bound
26//#include "CbcModel.hpp"
27//#include "IpCbcExtraData.hpp"
28
29//#include "CbcCompareUser.hpp"
30//#include "CbcCompareActual.hpp"
31//#include "CbcBranchUser.hpp"
32
33using namespace Ipopt;
34
35///Register options
36static void
37register_general_options
38(Ipopt::SmartPtr<Ipopt::RegisteredOptions> roptions)
39{
40  roptions->SetRegisteringCategory("bonmin output options");
41
42  roptions->AddBoundedIntegerOption("bb_log_level",
43      "specify main branch-and-bound log level.",
44      0,3,1,
45      "Set the level of output of the branch-and-bound : "
46      "0 - none, 1 - minimal, 2 - normal low, 3 - normal high"
47                                   );
48
49  roptions->AddLowerBoundedIntegerOption("bb_log_interval",
50      "Interval at which node level output is printed.",
51      0,100,
52      "Set the interval (in terms of number of nodes) at which "
53      "a log on node resolutions (consisting of lower and upper bounds) is given.");
54
55  roptions->AddBoundedIntegerOption("lp_log_level",
56      "specify LP log level.",
57      0,4,0,
58      "Set the level of output of the linear programming sub-solver in B-Hyb or B-QG : "
59      "0 - none, 1 - minimal, 2 - normal low, 3 - normal high, 4 - verbose"
60                                   );
61
62  roptions->AddBoundedIntegerOption("milp_log_level",
63      "specify MILP subsolver log level.",
64      0,3,0,
65      "Set the level of output of the MILP subsolver in OA : "
66      "0 - none, 1 - minimal, 2 - normal low, 3 - normal high"
67                                   );
68
69  roptions->AddBoundedIntegerOption("oa_log_level",
70      "specify OA iterations log level.",
71      0,2,1,
72      "Set the level of output of OA decomposition solver : "
73      "0 - none, 1 - normal, 2 - verbose"
74                                   );
75
76  roptions->AddLowerBoundedNumberOption("oa_log_frequency",
77      "display an update on lower and upper bounds in OA every n seconds",
78      0.,1.,100.,
79      "");
80
81  roptions->AddBoundedIntegerOption("nlp_log_level",
82      "specify NLP solver interface log level (independent from ipopt print_level).",
83      0,2,1,
84      "Set the level of output of the IpoptInterface : "
85      "0 - none, 1 - normal, 2 - verbose"
86                                   );
87
88  roptions->SetRegisteringCategory("bonmin branch-and-bound options");
89
90  roptions->AddStringOption4("algorithm",
91      "Choice of the algorithm.",
92      "B-Hyb",
93      "B-BB","simple branch-and-bound algorithm,",
94      "B-OA","OA Decomposition algorithm,",
95      "B-QG","Quesada and Grossmann branch-and-cut algorithm,",
96      "B-Hyb","hybrid outer approximation based branch-and-cut.",
97      "This will preset default values for most options of bonmin but depending on which algorithm "
98      "some of these can be changed.");
99  roptions->AddLowerBoundedNumberOption("time_limit",
100      "Set the global maximum computation time (in secs) for the algorithm.",
101      0.,1,1e10,
102      "");
103
104  roptions->AddLowerBoundedIntegerOption("node_limit",
105      "Set the maximum number of nodes explored in the branch-and-bound search.",
106      0,INT_MAX,
107      "");
108
109  roptions->AddBoundedNumberOption("integer_tolerance",
110      "Set integer tolerance.",
111      0.,1,.5,1,1e-06,
112      "Any number within that value of an integer is considered integer.");
113
114  roptions->AddBoundedNumberOption("allowable_gap",
115      "Specify the value of absolute gap under which the algorithm stops.",
116      -1.e20,0,1.e20,0,0.,
117      "Stop the tree search when the gap between the objective value of the best known solution"
118      " and the best bound on the objective of any solution is less than this.");
119
120  roptions->AddBoundedNumberOption("allowable_fraction_gap",
121      "Specify the value of relative gap under which the algorithm stops.",
122      -1.e20,0,1.e20,0,0.0,
123      "Stop the tree search when the gap between the objective value of the best known solution "
124      "and the best bound on the objective of any solution is less than this "
125      "fraction of the absolute value of the best known solution value.");
126
127  roptions->AddBoundedNumberOption("cutoff",
128      "Specify cutoff value.",
129      -1e100,0,1e100,0,1e100,
130      "cutoff should be the value of a feasible solution known by the user "
131      "(if any). The algorithm will only look for solutions better than cutoof.");
132
133
134  roptions->AddBoundedNumberOption("cutoff_decr",
135      "Specify cutoff decrement.",
136      -1e10,0,1e10,0,1e-05,
137      "Specify the amount by which cutoff is decremented below "
138      "a new best upper-bound"
139      " (usually a small postive value but in non-convex problems it may be a negative value).");
140
141  roptions->AddStringOption4("nodeselect_stra",
142      "Choose the node selection strategy.",
143      "best-bound",
144      "best-bound", "choose node whith the smallest bound,",
145      "depth-first", "Perform depth first search,",
146      "breadth-first", "Perform breadth first search,",
147      "dynamic", "Cbc dynamic strategy (starts with a depth first search and turn to best bound after 3 "
148      "integer feasible solutions have been found).",
149      "Choose the strategy for selecting the next node to be processed.");
150
151  roptions->AddLowerBoundedIntegerOption("number_strong_branch",
152      "Choose the maximum number of variables considered for strong branching.",
153      0,20,
154      "Set the number of variables on which to do strong branching.");
155
156  roptions->AddLowerBoundedIntegerOption
157  ("number_before_trust",
158   "Set the number of branches on a variable before its pseudo costs are to be believed "
159   "in dynamic strong branching.",
160   0,8,
161   "A value of 0 disables dynamic strong branching.");
162
163  roptions->AddStringOption3("warm_start",
164      "Select the warm start method",
165      "optimum",
166      "none","No warm start",
167      "optimum","Warm start with direct parent optimum",
168      "interior_point","Warm start with an interior point of direct parent",
169      "This will affect the function IpoptInterface::getWarmStart(), and as a consequence the wam starting in the various algorithms.");
170
171  roptions->AddStringOption2("sos_constraints",
172                             "Wether or not to activate SOS constraints.",
173                             "enable",
174                             "enable","",
175                             "disable","",
176                             "(only type 1 SOS are supported at the moment)");
177
178  roptions->SetRegisteringCategory("bonmin options for robustness");
179
180  roptions->AddLowerBoundedNumberOption("max_random_point_radius",
181      "Set max value r for coordinate of a random point.",
182      0.,1,1e5,
183      "When picking a random point coordinate i will be in the interval [min(max(l,-r),u-r), max(min(u,r),l+r)] "
184      "(where l is the lower bound for the variable and u is its upper bound)");
185
186  roptions->AddLowerBoundedIntegerOption
187  ("max_consecutive_failures",
188   "(temporarily removed) Number $n$ of consecutive unsolved problems before aborting a branch of the tree.",
189   0,10,
190   "When $n > 0$, continue exploring a branch of the tree until $n$ "
191   "consecutive problems in the branch are unsolved (we call unsolved a problem for which Ipopt can not "
192   "guarantee optimality within the specified tolerances).");
193
194  roptions->AddLowerBoundedIntegerOption
195  ("num_iterations_suspect",
196   "Number of iterations over which a node is considered \"suspect\" (for debugging purposes only, see detailed documentation).",
197   -1,-1,
198   "When the number of iterations to solve a node is above this number, the subproblem at this"
199   " node is considered to be suspect and it will be outputed in a file (set to -1 to deactivate this).");
200
201  roptions->AddStringOption2("nlp_failure_behavior",
202      "Set the behavior when an NLP or a series of NLP are unsolved by Ipopt (we call unsolved an NLP for which Ipopt is not "
203      "able to guarantee optimality within the specified tolerances).",
204      "stop",
205      "stop", "Stop when failure happens.",
206      "fathom", "Continue when failure happens.",
207      "If set to \"fathom\", the algorithm will fathom the node when Ipopt fails to find a solution to the nlp "
208      "at that node whithin the specified tolerances. "
209      "The algorithm then becomes a heuristic, and the user will be warned that the solution might not be optimal.");
210
211  roptions->AddLowerBoundedIntegerOption("num_retry_unsolved_random_point",
212      "Number $k$ of times that the algorithm will try to resolve an unsolved NLP with a random starting point "
213      "(we call unsolved an NLP for which Ipopt is not "
214      "able to guarantee optimality within the specified tolerances).",
215      0,0,
216      "When Ipopt fails to solve a continuous NLP sub-problem, if $k > 0$, the algorithm will "
217      "try again to solve the failed NLP with $k$ new randomly chosen starting points "
218      " or until the problem is solved with success.");
219
220
221  roptions->SetRegisteringCategory("bonmin options for non-convex problems");
222  roptions->AddLowerBoundedIntegerOption("max_consecutive_infeasible",
223      "Number of consecutive infeasible subproblems before aborting a"
224      " branch.",
225      0,0,
226      "Will continue exploring a branch of the tree until \"max_consecutive_infeasible\""
227      "consecutive problems are infeasibles by the NLP sub-solver.");
228
229  roptions->AddLowerBoundedIntegerOption("num_resolve_at_root",
230      "Number $k$ of tries to resolve the root node with different starting points.",
231      0,0,
232      "The algorithm will solve the root node with $k$ random starting points"
233      " and will keep the best local optimum found.");
234
235  roptions->AddLowerBoundedIntegerOption("num_resolve_at_node",
236      "Number $k$ of tries to resolve a node (other than the root) of the tree with different starting point.",
237      0,0,
238      "The algorithm will solve all the nodes with $k$ different random starting points "
239      "and will keep the best local optimum found.");
240
241
242}
243
244static void register_OA_options
245(Ipopt::SmartPtr<Ipopt::RegisteredOptions> roptions)
246{
247  roptions->SetRegisteringCategory("bonmin options : B-Hyb specific options");
248
249  roptions->AddLowerBoundedIntegerOption("nlp_solve_frequency",
250      "Specify the frequency (in terms of nodes) at which NLP relaxations are solved in B-Hyb.",
251      0,10,
252      "A frequency of 0 amounts to to never solve the NLP relaxation.");
253
254  roptions->AddLowerBoundedNumberOption("oa_dec_time_limit",
255      "Specify the maximum number of seconds spent overall in OA decomposition iterations.",
256      0.,0,30.,
257      "");
258  roptions->AddLowerBoundedNumberOption("tiny_element","Value for tiny element in OA cut",
259      -0.,0,1e-08,
260      "We will remove \"cleanly\" (by relaxing cut) an element lower"
261      " than this.");
262
263  roptions->AddLowerBoundedNumberOption("very_tiny_element","Value for very tiny element in OA cut",
264      -0.,0,1e-17,
265      "Algorithm will take the risk of neglecting an element lower"
266      " than this.");
267
268  roptions->AddLowerBoundedIntegerOption("Gomory_cuts",
269      "Frequency k (in terms of nodes) for generating Gomory cuts in branch-and-cut.",
270      -100,-5,
271      "If k > 0, cuts are generated every k nodes, if -99 < k < 0 cuts are generated every -k nodes but "
272      "Cbc may decide to stop generating cuts, if not enough are generated at the root node, "
273      "if k=-99 generate cuts only at the root node, if k=0 or 100 do not generate cuts.");
274  roptions->AddLowerBoundedIntegerOption("probing_cuts",
275      "Frequency (in terms of nodes) for generating probing cuts in branch-and-cut",
276      -100,-5,
277      "If k > 0, cuts are generated every k nodes, if -99 < k < 0 cuts are generated every -k nodes but "
278      "Cbc may decide to stop generating cuts, if not enough are generated at the root node, "
279      "if k=-99 generate cuts only at the root node, if k=0 or 100 do not generate cuts.");
280
281  roptions->AddLowerBoundedIntegerOption("cover_cuts",
282      "Frequency (in terms of nodes) for generating cover cuts in branch-and-cut",
283      -100,-5,
284      "If k > 0, cuts are generated every k nodes, if -99 < k < 0 cuts are generated every -k nodes but "
285      "Cbc may decide to stop generating cuts, if not enough are generated at the root node, "
286      "if k=-99 generate cuts only at the root node, if k=0 or 100 do not generate cuts.");
287
288  roptions->AddLowerBoundedIntegerOption("mir_cuts",
289      "Frequency (in terms of nodes) for generating MIR cuts in branch-and-cut",
290      -100,-5,
291      "If k > 0, cuts are generated every k nodes, if -99 < k < 0 cuts are generated every -k nodes but "
292      "Cbc may decide to stop generating cuts, if not enough are generated at the root node, "
293      "if k=-99 generate cuts only at the root node, if k=0 or 100 do not generate cuts.");
294
295}
296
297static void register_milp_sub_solver_options
298(Ipopt::SmartPtr<Ipopt::RegisteredOptions> roptions)
299{
300  roptions->SetRegisteringCategory("bonmin options : Options for milp subsolver in OA decomposition");
301  roptions->AddStringOption3("milp_subsolver",
302      "Choose the subsolver to solve MILPs sub-problems in OA decompositions.",
303      "Cbc_D",
304      "Cbc_D","Coin Branch and Cut with its default",
305      "Cbc_Par", "Coin Branch and Cut with passed parameters",
306      "Cplex","Ilog Cplex",
307      " To use Cplex, a valid license is required and you should have compiled OsiCpx in COIN-OR  (see Osi documentation).");
308}
309
310///Register options
311void
312IpoptInterface::register_ALL_options
313(Ipopt::SmartPtr<Ipopt::RegisteredOptions> roptions)
314{
315  register_general_options(roptions);
316  register_OA_options(roptions);
317  register_milp_sub_solver_options(roptions);
318}
319
320void
321IpoptInterface::set_ipopt_minlp_default(SmartPtr<OptionsList> Options)
322{
323  Options->SetNumericValue("gamma_phi", 1e-8, true, true);
324  Options->SetNumericValue("gamma_theta", 1e-4, true, true);
325  Options->SetNumericValue("required_infeasibility_reduction", 0.1, true, true);
326  Options->SetStringValue("expect_infeasible_problem","yes", true, true);
327  Options->SetStringValue("mu_strategy", "adaptive", true, true);
328  Options->SetStringValue("mu_oracle","probing", true, true);
329  Options->SetIntegerValue("print_level",1, true, true);
330}
331
332
333IpoptInterface::Messages::Messages
334():CoinMessages((int)IPOTPINTERFACE_DUMMY_END)
335{
336  strcpy(source_ ,"IpOp");
337  addMessage(SOLUTION_FOUND, CoinOneMessage
338      (1,2,"After %d tries found a solution of %g "
339       "(previous best %g)."));
340  addMessage(INFEASIBLE_SOLUTION_FOUND, CoinOneMessage
341      (2,2,"After %d tries found an solution of %g "
342       "infeasible problem."));
343
344  addMessage(UNSOLVED_PROBLEM_FOUND, CoinOneMessage
345      (3,2,"After %d tries found an solution of %g "
346       "unsolved problem."));
347  addMessage(WARN_SUCCESS_WS,CoinOneMessage
348      (3009,2,
349       "Problem not solved with warm start but "
350       "solved without"));
351
352  addMessage(WARNING_RESOLVING,CoinOneMessage
353      (3010,2,
354       "Trying to resolve NLP with different starting "
355       "point (%d attempts)."));
356  addMessage(WARN_SUCCESS_RANDOM, CoinOneMessage
357      (3011,1,
358       "Problem initially not solved but solved with "
359       "a random starting point (success on %d attempt)"));
360  addMessage(WARN_CONTINUING_ON_FAILURE, CoinOneMessage
361      (3017,1,
362       "Warning : continuing branching, while there are "
363       "unrecovered failures at nodes"));
364
365  addMessage(SUSPECT_PROBLEM, CoinOneMessage
366      (4,2,"NLP number %d is suspect (see bounds and start file)"));
367  addMessage(IPOPT_SUMMARY, CoinOneMessage
368      (6,2,"Ipopt return (for %s): status %2d, iter count %4d, time %g"));
369  addMessage(BETTER_SOL, CoinOneMessage
370      (7,2,"Solution of value %g found on %d'th attempt"));
371
372  addMessage(LOG_HEAD, CoinOneMessage
373      (8,1,"\n          "
374       "    Num      Status      Obj             It       time"));
375  addMessage(LOG_FIRST_LINE, CoinOneMessage
376      (9,1,
377       "    %-8d %-11s %-14g %-8d %-3g"));
378  addMessage(LOG_LINE, CoinOneMessage
379      (10,1,
380       " %c  r%-7d %-11s %-14g %-8d %-3g"));
381
382  addMessage(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, CoinOneMessage
383      (3012,1,"resolve called before any call to initialSolve"
384       " can not use warm starts."));
385
386}
387bool IpoptInterface::hasPrintedOptions=0;
388
389////////////////////////////////////////////////////////////////////
390// Constructors and desctructors                                  //
391////////////////////////////////////////////////////////////////////
392/// Default Constructor
393IpoptInterface::IpoptInterface():
394    OsiSolverInterface(),
395    tminlp_(NULL),
396    problem_(NULL),
397    app_(NULL),
398    rowsense_(NULL),
399    rhs_(NULL),
400    rowrange_(NULL),
401    reducedCosts_(NULL),
402    OsiDualObjectiveLimit_(1e200),
403//   optimization_status_(HasNotBeenOptimized),
404    varNames_(NULL),
405    hasVarNamesFile_(true),
406    nCallOptimizeTNLP_(0),
407    totalNlpSolveTime_(0),
408    totalIterations_(0),
409    maxRandomRadius_(1e08),
410    pushValue_(1e-02),
411    numRetryInitial_(-1),
412    numRetryResolve_(-1),
413    numRetryUnsolved_(1),
414    ipoptIMessages_(),
415    pretendFailIsInfeasible_(0),
416    hasContinuedAfterNlpFailure_(false),
417    numIterationSuspect_(-1),
418    hasBeenOptimized_(false),
419    warmStartStrategy_(1),
420    obj_(NULL),
421    feasibilityProblem_(NULL),
422    jRow_(NULL),
423    jCol_(NULL),
424    jValues_(NULL),
425    nnz_jac(0),
426    constTypes_(NULL),
427//    constTypesNum_(NULL),
428    nLinear_(0),
429    nNonLinear_(0),
430    tiny_(1e-8),
431    veryTiny_(1e-20),
432    infty_(1e100)
433{}
434
435
436/** Constructor with given IpSolver and TMINLP */
437IpoptInterface::IpoptInterface (Ipopt::SmartPtr<Ipopt::TMINLP> tminlp):
438    OsiSolverInterface(),
439    tminlp_(tminlp),
440    problem_(NULL),
441    app_(NULL),
442    rowsense_(NULL),
443    rhs_(NULL),
444    rowrange_(NULL),
445    reducedCosts_(NULL),
446    OsiDualObjectiveLimit_(1e200),
447//    optimization_status_(HasNotBeenOptimized),
448    varNames_(NULL),
449    hasVarNamesFile_(true),
450    nCallOptimizeTNLP_(0),
451    totalNlpSolveTime_(0),
452    totalIterations_(0),
453    maxRandomRadius_(1e08),
454    pushValue_(1e-02),
455    numRetryInitial_(-1),
456    numRetryResolve_(-1),
457    numRetryUnsolved_(1),
458    ipoptIMessages_(),
459    pretendFailIsInfeasible_(false),
460    hasContinuedAfterNlpFailure_(false),
461    numIterationSuspect_(-1),
462    hasBeenOptimized_(false),
463    warmStartStrategy_(1),
464    obj_(NULL),
465    feasibilityProblem_(NULL),
466    jRow_(NULL),
467    jCol_(NULL),
468    jValues_(NULL),
469    nnz_jac(0),
470    constTypes_(NULL),
471//    constTypesNum_(NULL),
472    nLinear_(0),
473    nNonLinear_(0),
474    tiny_(1e-08),
475    veryTiny_(1e-17),
476    infty_(1e100)
477{
478  assert(IsValid(tminlp));
479  app_ = new Ipopt::IpoptApplication();
480
481  SmartPtr<RegisteredOptions> roptions = app_->RegOptions();
482  register_ALL_options(roptions);
483  app_->Initialize("");
484  extractInterfaceParams();
485  set_ipopt_minlp_default(app_->Options());
486  // set the default options... expect_infeasible, etc...
487  problem_ = new Ipopt::TMINLP2TNLP(tminlp_, *app_->Options());
488  feasibilityProblem_ =
489    new Ipopt::TNLP2FPNLP
490    (Ipopt::SmartPtr<Ipopt::TNLP>(Ipopt::GetRawPtr(problem_)));
491
492}
493
494void
495IpoptInterface::readOptionFile(const char * fileName)
496{
497  app_->Initialize(fileName);
498  extractInterfaceParams();
499  set_ipopt_minlp_default(app_->Options());
500}
501void
502IpoptInterface::extractInterfaceParams()
503{
504  if (IsValid(app_)) {
505    app_->Options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,"bonmin.");
506    app_->Options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,"bonmin.");
507    app_->Options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,"bonmin.");
508    app_->Options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,"bonmin.");
509    app_->Options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,"bonmin.");
510    app_->Options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,"bonmin.");
511    app_->Options()->GetEnumValue("warm_start",warmStartStrategy_,"bonmin.");
512    app_->Options()->GetNumericValue
513    ("warm_start_bound_frac" ,pushValue_,"bonmin.");
514    app_->Options()->GetNumericValue("tiny_element",tiny_,"bonmin.");
515    app_->Options()->GetNumericValue("very_tiny_element",veryTiny_,"bonmin.");
516    double lower_infty, upper_infty;
517    app_->Options()->GetNumericValue("nlp_lower_bound_inf",lower_infty,"");
518    app_->Options()->GetNumericValue("nlp_upper_bound_inf",upper_infty,"");
519    infty_ = min(fabs(lower_infty), fabs(upper_infty));
520  }
521}
522
523/// Clone
524OsiSolverInterface * IpoptInterface::clone(bool CopyData) const
525{
526  //    debugMessage("IpoptInterface::clone(%d)\n", CopyData);
527  static int debug_no_clone = 0;
528  debug_no_clone++;
529  //    assert(debug_no_clone==1);
530  return new IpoptInterface(*this);
531}
532
533/// Copy constructor
534IpoptInterface::IpoptInterface (const IpoptInterface &source):
535    OsiSolverInterface(source),
536    tminlp_(source.tminlp_),
537    problem_(NULL),
538    app_(NULL),
539    rowsense_(NULL),
540    rhs_(NULL),
541    rowrange_(NULL),
542    reducedCosts_(NULL),
543    OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
544    optimization_status_(source.optimization_status_),
545    varNames_(NULL),
546    hasVarNamesFile_(source.hasVarNamesFile_),
547    nCallOptimizeTNLP_(0),
548    totalNlpSolveTime_(0),
549    totalIterations_(0),
550    maxRandomRadius_(source.maxRandomRadius_),
551    pushValue_(source.pushValue_),
552    numRetryInitial_(source.numRetryInitial_),
553    numRetryResolve_(source.numRetryResolve_),
554    numRetryUnsolved_(source.numRetryUnsolved_),
555    ipoptIMessages_(),
556    pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
557    hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
558    numIterationSuspect_(source.numIterationSuspect_),
559    hasBeenOptimized_(source.hasBeenOptimized_),
560    warmStartStrategy_(source.warmStartStrategy_),
561    obj_(NULL),
562    feasibilityProblem_(NULL),
563    jRow_(NULL),
564    jCol_(NULL),
565    jValues_(NULL),
566    nnz_jac(source.nnz_jac),
567    constTypes_(NULL),
568//    constTypesNum_(NULL),
569    nLinear_(0),
570    nNonLinear_(0),
571    tiny_(source.tiny_),
572    veryTiny_(source.veryTiny_),
573    infty_(source.infty_)
574{
575  // Copy options from old application
576  if(IsValid(source.tminlp_)) {
577    app_ = new Ipopt::IpoptApplication();
578
579    SmartPtr<RegisteredOptions> roptions = app_->RegOptions();
580    register_ALL_options(roptions);
581    // Copy the options
582    *app_->Options()=*source.app_->Options();
583
584    //    extractInterfaceParams();
585
586    problem_ = new Ipopt::TMINLP2TNLP(tminlp_, *app_->Options());
587    problem_->copyUserModification(*source.problem_);
588    pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
589  }
590  else {
591    throw SimpleError("Don't know how to copy an empty IpoptInterface.",
592        "copy constructor");
593  }
594
595  if(source.obj_) {
596    obj_ = new double[source.getNumCols()];
597    CoinCopyN(source.obj_, source.getNumCols(), obj_);
598  }
599  if(IsValid(source.tminlp_))
600    feasibilityProblem_ = new Ipopt::TNLP2FPNLP
601        (Ipopt::SmartPtr<Ipopt::TNLP>(Ipopt::GetRawPtr(problem_)));
602  else
603    throw CoinError("Don't know how to copy an empty IpoptOAInterface.",
604        "copy constructor","IpoptOAInterface");
605
606
607  if(source.jValues_!=NULL && source.jRow_ != NULL && source.jCol_ != NULL && nnz_jac>0) {
608    jValues_ = new double [nnz_jac];
609    jCol_    = new Ipopt::Index [nnz_jac];
610    jRow_    = new Ipopt::Index [nnz_jac];
611    CoinCopyN(source.jValues_ , nnz_jac,jValues_ );
612    CoinCopyN(source.jCol_    , nnz_jac,jCol_    );
613    CoinCopyN(source.jRow_    , nnz_jac,jRow_    );
614
615    if(source.constTypes_ != NULL) {
616      constTypes_ = new Ipopt::TMINLP::ConstraintType[getNumRows()];
617      CoinCopyN(source.constTypes_, getNumRows(), constTypes_);
618    }
619//    if(source.constTypesNum_ != NULL) {
620//      constTypesNum_ = new int[getNumRows()];
621//      CoinCopyN(source.constTypesNum_, getNumRows(), constTypesNum_);
622//    }
623  }
624  else if(nnz_jac > 0) {
625    throw CoinError("Arrays for storing jacobian are inconsistant.",
626        "copy constructor","IpoptOAInterface");
627  }
628  // Process the output options
629  app_->Initialize("");
630}
631
632
633Ipopt::SmartPtr<Ipopt::OptionsList> IpoptInterface::retrieve_options()
634{
635  if(!IsValid(app_)) {
636    std::cout<<"Can not parse options when no IpApplication has been created"<<std::endl;
637    return NULL;
638  }
639  else
640    return app_->Options();
641}
642
643/// Assignment operator
644IpoptInterface & IpoptInterface::operator=(const IpoptInterface& rhs)
645{
646  if(this!= &rhs) {
647    OsiSolverInterface::operator=(rhs);
648    OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
649    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
650    totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
651    totalIterations_ = rhs.totalIterations_;
652    maxRandomRadius_ = rhs.maxRandomRadius_;
653    hasVarNamesFile_ = rhs.hasVarNamesFile_;
654    pushValue_ = rhs.pushValue_;
655    optimization_status_ = rhs.optimization_status_;
656
657    if(IsValid(rhs.tminlp_)) {
658
659      tminlp_ = rhs.tminlp_;
660      app_ = new Ipopt::IpoptApplication();
661      problem_ = new Ipopt::TMINLP2TNLP(tminlp_, *app_->Options());
662
663      feasibilityProblem_ = new Ipopt::TNLP2FPNLP
664          (Ipopt::SmartPtr<Ipopt::TNLP>(Ipopt::GetRawPtr(problem_)));
665      nnz_jac = rhs.nnz_jac;
666
667      if(constTypes_ != NULL) {
668        delete [] constTypes_;
669        constTypes_ = NULL;
670      }
671      if(rhs.constTypes_ != NULL) {
672        constTypes_ = new Ipopt::TMINLP::ConstraintType[getNumRows()];
673        CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
674      }
675/*
676      if(constTypesNum_ != NULL) {
677        delete [] constTypesNum_;
678        constTypesNum_ = NULL;
679      }
680      if(rhs.constTypesNum_ != NULL) {
681        constTypesNum_ = new int[getNumRows()];
682        CoinCopyN(rhs.constTypesNum_, getNumRows(), constTypesNum_);
683      }
684*/
685      if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
686        jValues_ = new double [nnz_jac];
687        jCol_    = new Ipopt::Index [nnz_jac];
688        jRow_    = new Ipopt::Index [nnz_jac];
689        CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
690        CoinCopyN(rhs.jCol_    , nnz_jac,jCol_    );
691        CoinCopyN(rhs.jRow_    , nnz_jac,jRow_    );
692      }
693      else if(nnz_jac > 0) {
694        throw CoinError("Arrays for storing jacobian are inconsistant.",
695            "copy constructor",
696            "IpoptOAInterface");
697      }
698      tiny_ = rhs.tiny_;
699      veryTiny_ = rhs.veryTiny_;
700      infty_ = rhs.infty_;
701
702
703    }
704    else {
705      tminlp_ =NULL;
706      app_ = NULL;
707      problem_ = NULL;
708      feasibilityProblem_ = NULL;
709    }
710
711
712    if(obj_) {
713      delete [] obj_;
714      obj_ = NULL;
715    }
716    if(rhs.obj_) {
717      obj_ = new double[rhs.getNumCols()];
718      CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
719    }
720
721    delete [] varNames_;
722    varNames_ = NULL;
723
724    if(rhs.varNames_) {
725      rhs.varNames_ = new std::string[getNumCols()];
726      CoinCopyN(rhs.varNames_, getNumCols(), varNames_);
727    }
728
729    hasVarNamesFile_ = rhs.hasVarNamesFile_;
730
731    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
732    totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
733    totalIterations_ = rhs.totalIterations_;
734    maxRandomRadius_ = rhs.maxRandomRadius_;
735    pushValue_ = rhs.pushValue_;
736    numRetryInitial_ = rhs.numRetryInitial_;
737    numRetryResolve_ = rhs.numRetryResolve_;
738    numRetryUnsolved_ = rhs.numRetryUnsolved_;
739    pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
740    numIterationSuspect_ = rhs.numIterationSuspect_;
741
742    hasBeenOptimized_ = rhs.hasBeenOptimized_;
743
744    freeCachedData();
745  }
746  return *this;
747}
748
749/// Destructor
750IpoptInterface::~IpoptInterface ()
751{
752  freeCachedData();
753  delete [] jRow_;
754  delete [] jCol_;
755  delete [] jValues_;
756  delete [] constTypes_;
757//  delete [] constTypesNum_;
758}
759
760void
761IpoptInterface::freeCachedColRim()
762{
763  if(varNames_!=NULL) {
764    delete [] varNames_;
765    varNames_ = NULL;
766  }
767  if(reducedCosts_!=NULL) {
768    delete []  reducedCosts_;
769    reducedCosts_ = NULL;
770  }
771
772}
773
774void
775IpoptInterface::freeCachedRowRim()
776{
777  if(rowsense_!=NULL) {
778    delete []  rowsense_;
779    rowsense_ = NULL;
780  }
781  if(rhs_!=NULL) {
782    delete []  rhs_;
783    rhs_ = NULL;
784  }
785  if(rowrange_!=NULL) {
786    delete []  rowrange_;
787    rowrange_ = NULL;
788  }
789  //   if(dualsol_!=NULL)
790  //     {
791  //       delete []  dualsol_; dualsol_ = NULL;
792  //     }
793}
794
795void
796IpoptInterface::freeCachedData()
797{
798  freeCachedColRim();
799  freeCachedRowRim();
800}
801
802///////////////////////////////////////////////////////////////////
803// WarmStart Information                                                                           //
804///////////////////////////////////////////////////////////////////
805
806/// Get warmstarting information
807CoinWarmStart*
808IpoptInterface::getWarmStart() const
809{
810  if(warmStartStrategy_) {
811    if(warmStartStrategy_==2) {
812      SmartPtr<IpoptInteriorWarmStarter> warm_starter =
813        SmartPtr<IpoptInteriorWarmStarter>(problem_->GetWarmStarter());
814      return new IpoptWarmStart(*this, warm_starter);
815    }
816    else  return new IpoptWarmStart(*this, NULL);
817  }
818  else
819    return new IpoptWarmStart(getNumCols(), getNumRows());
820}
821
822
823bool
824IpoptInterface::setWarmStart(const CoinWarmStart* warmstart)
825{
826  if(!warmstart || !warmStartStrategy_)
827    return 0;
828  hasBeenOptimized_ = false;
829  const IpoptWarmStart * ws = dynamic_cast<const IpoptWarmStart*> (warmstart);
830  if(ws->empty())//reset initial point and leave
831  {
832    unsetWarmStartOptions();
833    return 1;
834  }
835  setWarmStartOptions();
836  int numcols = getNumCols();
837  int numrows = getNumRows();
838  const double * colLow = getColLower();
839  const double * colUp = getColUpper();
840  for(int i = 0; i < numcols ; i++) {
841    CoinWarmStartBasis::Status status = ws->getStructStatus(i);
842    if(status == CoinWarmStartBasis::atLowerBound) {
843      problem_->setxInit(i,colLow[i]);
844      problem_->setDualInit(i + numcols + numrows,0.);
845    }
846    else if(status == CoinWarmStartBasis::atUpperBound) {
847      problem_->setxInit(i,colUp[i]);
848      problem_->setDualInit(i + numrows,0.);
849    }
850    else {
851      problem_->setDualInit(i + numrows,0.);
852      problem_->setDualInit(i + numcols + numrows, 0.);
853    }
854  }
855  for(int i = 0; i < numrows ; i++) {
856    CoinWarmStartBasis::Status status = ws->getArtifStatus(i);
857    if(status == CoinWarmStartBasis::atLowerBound) {
858      problem_->setDualInit(i,0.);
859    }
860  }
861  int nElem = ws->values()->getNumElements();
862  const int * inds = ws->values()->getIndices();
863  const double * elems = ws->values()->getElements();
864
865  for(int i = 0 ; i < nElem ; i++) {
866    problem_->setxInit(inds[i],elems[i]);
867  }
868
869  if(IsValid(ws->warm_starter()))
870    problem_->SetWarmStarter(ws->warm_starter());
871  return 1;
872}
873
874static const char * OPT_SYMB="OPT";
875static const char * FAILED_SYMB="FAILED";
876static const char * INFEAS_SYMB="INFEAS";
877
878void
879IpoptInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
880    const char * whereFrom)
881{
882  totalNlpSolveTime_-=CoinCpuTime();
883  if(warmStarted)
884    optimization_status_ = app_->ReOptimizeTNLP(GetRawPtr(problem_));
885  else
886    optimization_status_ = app_->OptimizeTNLP(GetRawPtr(problem_));
887  totalNlpSolveTime_+=CoinCpuTime();
888  nCallOptimizeTNLP_++;
889  hasBeenOptimized_ = true;
890
891  //Options should have been printed if not done already turn off Ipopt output
892  if(!hasPrintedOptions) {
893    hasPrintedOptions = 1;
894    //app_->Options()->SetIntegerValue("print_level",0, true, true);
895    app_->Options()->SetStringValue("print_user_options","no", false, true);
896  }
897
898  const SmartPtr<SolveStatistics>  stats = app_->Statistics();
899
900#if 1
901  if(optimization_status_ == -10)//Too few degrees of freedom
902    {
903      std::cout<<"Too few degrees of freedom...."<<std::endl;
904      int numrows = getNumRows();
905      int numcols = getNumCols();
906
907      const double * colLower = getColLower();
908      const double * colUpper = getColUpper();
909
910
911      const double * rowLower = getRowLower();
912      const double * rowUpper = getRowUpper();
913
914      int numberFixed = 0;
915      for(int i = 0 ; i < numcols ; i++)
916        {
917          if(colUpper[i] - colLower[i] <= INT_BIAS)
918            {
919              numberFixed++;
920            }
921        }
922      int numberEqualities = 0;
923      for(int i = 0 ; i < numrows ; i++)
924        {
925          if(rowUpper[i] - rowLower[i] <= INT_BIAS)
926            {
927              numberEqualities++;
928            }     
929        }
930      if(numcols - numberFixed > numberEqualities)
931        {
932          throw UnsolvedError(optimization_status_);
933        }
934      double * saveColLow = CoinCopyOfArray(getColLower(), getNumCols());
935      double * saveColUp = CoinCopyOfArray(getColUpper(), getNumCols());
936
937      for(int i = 0; i < numcols && numcols - numberFixed <= numberEqualities ; i++)
938        {
939          if(colUpper[i] - colLower[i] <= INT_BIAS)
940            {
941              setColLower(i, saveColLow[i]-1e-06);
942              setColUpper(i, saveColLow[i]+1e-06);
943              numberFixed--;
944            }
945        }
946      solveAndCheckErrors(warmStarted, throwOnFailure, whereFrom);
947      //restore
948      for(int i = 0; i < numcols && numcols - numberFixed < numrows ; i++)
949        {
950          problem_->SetVariableLowerBound(i,saveColLow[i]);
951          problem_->SetVariableUpperBound(i,saveColUp[i]);
952        }
953      delete [] saveColLow;
954      delete [] saveColUp;
955      return;
956    }
957  else 
958#endif
959    if(optimization_status_ < -9)//Ipopt failed and the error can not be recovered, throw it
960  {
961    throw UnsolvedError(optimization_status_);
962  }
963
964  if(IsValid(stats)) {
965    totalIterations_ += stats->IterationCount();
966  }
967  else if (throwOnFailure)//something failed throw
968  {
969    throw SimpleError("No statistics available from Ipopt",whereFrom);
970  }
971  else return;
972
973  messageHandler()->message(IPOPT_SUMMARY, ipoptIMessages_)
974  <<whereFrom<<optimization_status_<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
975
976  if((nCallOptimizeTNLP_ % 20) == 1)
977    messageHandler()->message(LOG_HEAD, ipoptIMessages_)<<CoinMessageEol;
978
979
980  if (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) {
981    messageHandler()->message(SUSPECT_PROBLEM,
982        ipoptIMessages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
983    std::string subProbName;
984    getStrParam(OsiProbName, subProbName);
985    std::ostringstream os;
986    os<<"_"<<nCallOptimizeTNLP_;
987    subProbName+=os.str();
988    problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
989  }
990
991}
992
993////////////////////////////////////////////////////////////////////
994// Solve Methods                                                  //
995////////////////////////////////////////////////////////////////////
996/// Solve initial continuous relaxation
997void IpoptInterface::initialSolve()
998{
999  assert(IsValid(app_));
1000  assert(IsValid(problem_));
1001  if(!problem_->checkZeroDimension(optimization_status_)) {
1002
1003    if(!hasPrintedOptions) {
1004      int printOptions;
1005      app_->Options()->GetEnumValue("print_user_options",printOptions,"bonmin.");
1006      if(printOptions)
1007        app_->Options()->SetStringValue("print_user_options","yes");
1008    }
1009    solveAndCheckErrors(0,1,"initialSolve");
1010
1011    //Options should have been printed if not done already turn off Ipopt output
1012    if(!hasPrintedOptions) {
1013      hasPrintedOptions = 1;
1014      app_->Options()->SetStringValue("print_user_options","no");
1015      app_->Options()->SetIntegerValue("print_level",0);
1016    }
1017
1018    const char * status=OPT_SYMB;
1019    ;
1020    if(isAbandoned()) status=FAILED_SYMB;
1021    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1022    messageHandler()->message(LOG_FIRST_LINE, ipoptIMessages_)<<nCallOptimizeTNLP_
1023    <<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1024
1025    if(isAbandoned()) {
1026      resolveForRobustness(numRetryUnsolved_);
1027    }
1028    else if(numRetryInitial_)
1029    {
1030      resolveForCost(numRetryInitial_);
1031      /** Only do it once at the root.*/
1032      numRetryInitial_ = 0;
1033    }
1034  }
1035  else
1036    hasBeenOptimized_ = true;
1037}
1038
1039/** Resolve the continuous relaxation after problem modification.
1040 * \note for Ipopt, same as resolve */
1041void
1042IpoptInterface::resolve()
1043{
1044  if (!IsValid(app_->Statistics())) {
1045    messageHandler()->message(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, ipoptIMessages_)
1046    <<CoinMessageEol;
1047    initialSolve();
1048    return;
1049  }
1050  assert(IsValid(app_));
1051  assert(IsValid(problem_));
1052  if (INT_BIAS > 0.) {
1053    app_->Options()->SetStringValue("warm_start_same_structure", "yes");
1054  }
1055  else {
1056    app_->Options()->SetStringValue("warm_start_same_structure", "no");
1057  }
1058
1059  setWarmStartOptions();
1060
1061  if(!problem_->checkZeroDimension(optimization_status_)) {
1062    solveAndCheckErrors(1,1,"resolve");
1063
1064    const char * status=OPT_SYMB;
1065    ;
1066    if(isAbandoned()) status=FAILED_SYMB;
1067    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1068    messageHandler()->message(LOG_FIRST_LINE, ipoptIMessages_)<<nCallOptimizeTNLP_
1069    <<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1070
1071    if(isAbandoned()) {
1072      resolveForRobustness(numRetryUnsolved_);
1073    }
1074    else if(numRetryResolve_)
1075      resolveForCost(numRetryResolve_);
1076  }
1077  else
1078    hasBeenOptimized_ = true;
1079}
1080
1081void
1082IpoptInterface::resolveForCost(int numsolve)
1083{
1084
1085  /** Save current optimum. */
1086  double * point = new double[getNumCols()*3+ getNumRows()];
1087  double bestBound = getObjValue();
1088  CoinCopyN(getColSolution(),
1089      getNumCols(), point);
1090  CoinCopyN(getRowPrice(),
1091      2*getNumCols()+ getNumRows(),
1092      &point[getNumCols()]);
1093
1094  if(isProvenOptimal())
1095    messageHandler()->message(SOLUTION_FOUND,
1096        ipoptIMessages_)
1097    <<1<<getObjValue()<<bestBound
1098    <<CoinMessageEol;
1099  else
1100    messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
1101        ipoptIMessages_)
1102    <<1
1103    <<CoinMessageEol;
1104  for(int f = 0; f < numsolve ; f++) {
1105    messageHandler()->message(WARNING_RESOLVING,
1106        ipoptIMessages_)
1107    <<f+1<< CoinMessageEol ;
1108    randomStartingPoint();
1109    solveAndCheckErrors(0,0,"resolveForCost");
1110
1111
1112    const SmartPtr<SolveStatistics>  stats = app_->Statistics();
1113    if(IsValid(stats)) {
1114      totalIterations_ += stats->IterationCount();
1115    }
1116    else//something failed (random point generation is not that clever and this can happen (for ex asked to compute sqrt(-1))
1117      //just exit this procedure
1118    {
1119      return;
1120      //throw SimpleError("No statistics available from Ipopt","resolveForCost");
1121    }
1122
1123
1124
1125
1126
1127    const char * status=OPT_SYMB;
1128    ;
1129    char c=' ';
1130    if(isAbandoned()) {
1131      status=FAILED_SYMB;
1132    }
1133    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1134
1135
1136    //Is solution better than previous
1137    if(isProvenOptimal() &&
1138        getObjValue()<bestBound) {
1139      c='*';
1140      messageHandler()->message(BETTER_SOL, ipoptIMessages_)<<getObjValue()<<f+1<< CoinMessageEol;
1141      CoinCopyN(getColSolution(),
1142          getNumCols(), point);
1143      CoinCopyN(getRowPrice(),
1144          2*getNumCols()+ getNumRows(),
1145          &point[getNumCols()]);
1146      bestBound = getObjValue();
1147    }
1148
1149    messageHandler()->message(LOG_LINE, ipoptIMessages_)
1150    <<c<<f+1<<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1151
1152
1153    if(isProvenOptimal())
1154      messageHandler()->message(SOLUTION_FOUND,
1155          ipoptIMessages_)
1156      <<f+2<<getObjValue()<<bestBound
1157      <<CoinMessageEol;
1158    else if(!isAbandoned())
1159      messageHandler()->message(UNSOLVED_PROBLEM_FOUND,
1160          ipoptIMessages_)
1161      <<f+2
1162      <<CoinMessageEol;
1163    else
1164      messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
1165          ipoptIMessages_)
1166      <<f+2
1167      <<CoinMessageEol;
1168  }
1169  setColSolution(point);
1170  setRowPrice(&point[getNumCols()]);
1171  setWarmStartOptions();
1172  delete [] point;
1173
1174  optimization_status_ = app_->ReOptimizeTNLP(GetRawPtr(problem_));
1175  hasBeenOptimized_ = true;
1176}
1177
1178void
1179IpoptInterface::resolveForRobustness(int numsolve)
1180{
1181  //std::cerr<<"Resolving the problem for robustness"<<std::endl;
1182  //First remove warm start point and resolve
1183  unsetWarmStartOptions();
1184  messageHandler()->message(WARNING_RESOLVING,
1185      ipoptIMessages_)
1186  <<1<< CoinMessageEol ;
1187  solveAndCheckErrors(0,0,"resolveForRobustness");
1188
1189
1190  const char * status=OPT_SYMB;
1191  ;
1192  char c='*';
1193  if(isAbandoned()) {
1194    status=FAILED_SYMB;
1195    c=' ';
1196  }
1197  else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1198  messageHandler()->message(LOG_LINE, ipoptIMessages_)
1199  <<c<<1<<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1200
1201
1202  if(!isAbandoned()) {
1203    messageHandler()->message(WARN_SUCCESS_WS,
1204        ipoptIMessages_)
1205    << CoinMessageEol ;
1206    return; //we won go on
1207  }
1208
1209  //still unsolved try again with different random starting points
1210  for(int f = 0; f < numsolve ; f++) {
1211    messageHandler()->message(WARNING_RESOLVING,
1212        ipoptIMessages_)
1213    <<f+2<< CoinMessageEol ;
1214
1215    randomStartingPoint();
1216    solveAndCheckErrors(0,0,"resolveForRobustness");
1217
1218
1219    messageHandler()->message(IPOPT_SUMMARY, ipoptIMessages_)
1220    <<"resolveForRobustness"<<optimization_status_<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1221
1222
1223    const char * status=OPT_SYMB;
1224    ;
1225    char c='*';
1226    if(isAbandoned()) {
1227      status=FAILED_SYMB;
1228      c=' ';
1229    }
1230    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1231    messageHandler()->message(LOG_LINE, ipoptIMessages_)
1232    <<c<<f+2<<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1233
1234
1235    if(!isAbandoned()) {
1236      messageHandler()->message(WARN_SUCCESS_RANDOM,
1237          ipoptIMessages_)
1238      <<f+2
1239      << CoinMessageEol ;
1240      return; //we have found a solution and continue
1241    }
1242  }
1243  if(pretendFailIsInfeasible_) {
1244    if(pretendFailIsInfeasible_ == 1) {
1245      messageHandler()->message(WARN_CONTINUING_ON_FAILURE,
1246          ipoptIMessages_)
1247      <<CoinMessageEol;
1248      hasContinuedAfterNlpFailure_ = 1;
1249    }
1250    return;
1251  }
1252  else {
1253    throw UnsolvedError(optimization_status_);
1254  }
1255}
1256
1257////////////////////////////////////////////////////////////////////
1258// Methods returning info on how the solution process terminated  //
1259////////////////////////////////////////////////////////////////////
1260/// Are there a numerical difficulties?
1261bool IpoptInterface::isAbandoned() const
1262{
1263  return (
1264        (optimization_status_==Ipopt::Maximum_Iterations_Exceeded)||
1265        (optimization_status_==Ipopt::Restoration_Failed)||
1266        (optimization_status_==Ipopt::Error_In_Step_Computation)||
1267        (optimization_status_==Ipopt::Not_Enough_Degrees_Of_Freedom)||
1268        (optimization_status_==Ipopt::Invalid_Problem_Definition)||
1269        (optimization_status_==Ipopt::Invalid_Option)||
1270        (optimization_status_==Ipopt::Invalid_Number_Detected)||
1271        (optimization_status_==Ipopt::Unrecoverable_Exception)||
1272        (optimization_status_==Ipopt::NonIpopt_Exception_Thrown)||
1273        (optimization_status_==Ipopt::Insufficient_Memory)||
1274        (optimization_status_==Ipopt::Internal_Error)
1275      );
1276}
1277
1278/// Is optimality proven?
1279bool IpoptInterface::isProvenOptimal() const
1280{
1281  if (optimization_status_==Ipopt::Search_Direction_Becomes_Too_Small) {
1282    std::cerr<<"Warning : need to verify that Search_Direction_Becomes_Too_Small is indeed Ok"<<std::endl;
1283  }
1284  return (optimization_status_==Ipopt::Solve_Succeeded ||
1285      optimization_status_==Ipopt::Search_Direction_Becomes_Too_Small ||
1286      optimization_status_==Ipopt::Solved_To_Acceptable_Level);
1287}
1288/// Is primal infeasiblity proven?
1289bool IpoptInterface::isProvenPrimalInfeasible() const
1290{
1291  return (optimization_status_ == Ipopt::Infeasible_Problem_Detected);
1292}
1293/// Is dual infeasiblity proven?
1294bool IpoptInterface::isProvenDualInfeasible() const
1295{
1296  throw SimpleError("Don't have this optimization status yet.",
1297      "isProvenDualInfeasible");
1298}
1299/// Is the given primal objective limit reached?
1300bool IpoptInterface::isPrimalObjectiveLimitReached() const
1301{
1302  std::cerr<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<std::endl;
1303  return 0;
1304}
1305/// Is the given dual objective limit reached?
1306bool IpoptInterface::isDualObjectiveLimitReached() const
1307{
1308  //  std::cerr<<"Warning : isDualObjectiveLimitReached not implemented yet"<<std::endl;
1309  return (optimization_status_==Ipopt::Diverging_Iterates);
1310
1311}
1312/// Iteration limit reached?
1313bool IpoptInterface::isIterationLimitReached() const
1314{
1315  return (optimization_status_==Ipopt::Maximum_Iterations_Exceeded);
1316  //  return (problem_->optimization_status()==Ipopt::Maximum_Iterations_Exceeded);
1317}
1318////////////////////////////////////////////////////////////////////
1319// Problem information methods                                    //
1320////////////////////////////////////////////////////////////////////
1321/// Get number of columns
1322int IpoptInterface::getNumCols() const
1323{
1324
1325  return problem_->num_variables();
1326}
1327
1328
1329/// Get number of rows
1330int
1331IpoptInterface::getNumRows() const
1332{
1333  return problem_->num_constraints();
1334}
1335
1336const double *
1337IpoptInterface::getColLower() const
1338{
1339  return problem_->x_l();
1340}
1341
1342const double *
1343IpoptInterface::getColUpper() const
1344{
1345  return problem_->x_u();
1346}
1347
1348void
1349IpoptInterface::readVarNames() const
1350{
1351  delete []varNames_;
1352  varNames_ = NULL;
1353  std::string probName;
1354  getStrParam(OsiProbName, probName);
1355  IpCbcColReader colRead(probName);
1356  if(colRead.readFile()) {
1357    varNames_ = new std::string[problem_->num_variables()];
1358    colRead.copyNames(varNames_,problem_->num_variables());
1359    hasVarNamesFile_ = true;
1360  }
1361  else
1362    hasVarNamesFile_ = false;
1363}
1364
1365///get name of a variable
1366const std::string *
1367IpoptInterface::getVarNames() const
1368{
1369  if(varNames_ == NULL && hasVarNamesFile_ ) {
1370    readVarNames();
1371  }
1372  return varNames_;
1373}
1374
1375void IpoptInterface::extractSenseRhsAndRange() const
1376{
1377  assert(rowsense_==NULL&&rhs_==NULL&&rowrange_==NULL);
1378  int numrows = problem_->num_constraints();
1379  const double * rowLower = getRowLower();
1380  const double * rowUpper = getRowUpper();
1381  rowsense_ = new char [numrows];
1382  rhs_ = new double [numrows];
1383  rowrange_ = new double [numrows];
1384  for(int i = 0 ; i < numrows ; i++) {
1385    rowrange_[i]=0.;
1386    convertBoundToSense(rowLower[i], rowUpper[i], rowsense_[i], rhs_[i], rowrange_[i]);
1387  }
1388}
1389
1390/** Get pointer to array[getNumRows()] of row constraint senses.
1391    <ul>
1392    <li>'L': <= constraint
1393    <li>'E': =  constraint
1394    <li>'G': >= constraint
1395    <li>'R': ranged constraint
1396    <li>'N': free constraint
1397    </ul>
1398*/
1399const char *
1400IpoptInterface::getRowSense() const
1401{
1402  if(rowsense_==NULL) {
1403    extractSenseRhsAndRange();
1404  }
1405  return rowsense_;
1406}
1407
1408/** Get pointer to array[getNumRows()] of rows right-hand sides
1409    <ul>
1410    <li> if rowsense()[i] == 'L' then rhs()[i] == rowupper()[i]
1411    <li> if rowsense()[i] == 'G' then rhs()[i] == rowlower()[i]
1412    <li> if rowsense()[i] == 'R' then rhs()[i] == rowupper()[i]
1413    <li> if rowsense()[i] == 'N' then rhs()[i] == 0.0
1414    </ul>
1415*/
1416const double *
1417IpoptInterface::getRightHandSide() const
1418{
1419  if(rhs_==NULL) {
1420    extractSenseRhsAndRange();
1421  }
1422  return rhs_;
1423}
1424
1425/** Get pointer to array[getNumRows()] of row ranges.
1426    <ul>
1427    <li> if rowsense()[i] == 'R' then
1428    rowrange()[i] == rowupper()[i] - rowlower()[i]
1429    <li> if rowsense()[i] != 'R' then
1430    rowrange()[i] is 0.0
1431    </ul>
1432*/
1433const double *
1434IpoptInterface::getRowRange() const
1435{
1436  if(rowrange_==NULL) {
1437    extractSenseRhsAndRange();
1438  }
1439  return rowrange_;
1440}
1441
1442const double *
1443IpoptInterface::getRowLower() const
1444{
1445  return problem_->g_l();
1446}
1447
1448const double *
1449IpoptInterface::getRowUpper() const
1450{
1451  return problem_->g_u();
1452}
1453
1454/// Return true if column is continuous
1455bool
1456IpoptInterface::isContinuous(int colNumber) const
1457{
1458  return (problem_->var_types()[colNumber]==Ipopt::TMINLP::CONTINUOUS);
1459}
1460
1461/// Return true if column is binary
1462bool
1463IpoptInterface::isBinary(int colNumber) const
1464{
1465  return (problem_->var_types()[colNumber]==Ipopt::TMINLP::BINARY);
1466}
1467
1468/** Return true if column is integer.
1469    Note: This function returns true if the the column
1470    is binary or a general integer.
1471*/
1472bool
1473IpoptInterface::isInteger(int colNumber) const
1474{
1475  return ((problem_->var_types()[colNumber]==Ipopt::TMINLP::BINARY)||
1476      (problem_->var_types()[colNumber]==Ipopt::TMINLP::INTEGER));
1477}
1478
1479/// Return true if column is general integer
1480bool
1481IpoptInterface::isIntegerNonBinary(int colNumber) const
1482{
1483  return (problem_->var_types()[colNumber]==Ipopt::TMINLP::INTEGER);
1484}
1485/// Return true if column is binary and not fixed at either bound
1486bool
1487IpoptInterface::isFreeBinary(int colNumber) const
1488{
1489  return ((problem_->var_types()[colNumber]==Ipopt::TMINLP::BINARY)
1490      &&((getColUpper()[colNumber]-getColLower()[colNumber]) > 1 - 1e-09));
1491}
1492
1493/// Get solver's value for infinity
1494double
1495IpoptInterface::getInfinity() const
1496{
1497  return infty_;
1498}
1499
1500/// Get pointer to array[getNumCols()] of primal solution vector
1501const double *
1502IpoptInterface::getColSolution() const
1503{
1504  if(hasBeenOptimized_)
1505    return problem_->x_sol();
1506  else
1507    return problem_->x_init();
1508}
1509
1510/// Get pointer to array[getNumRows()] of dual prices
1511const double *
1512IpoptInterface::getRowPrice() const
1513{
1514  if(hasBeenOptimized_)
1515    return problem_->duals_sol();
1516  else
1517    return problem_->duals_init();
1518}
1519
1520/// Get a pointer to array[getNumCols()] of reduced costs
1521const double *
1522IpoptInterface::getReducedCost() const
1523{
1524  std::cerr<<"WARNING : trying to access reduced cost in Ipopt always retrun 0"<<std::endl;
1525  if(reducedCosts_==NULL) {
1526    reducedCosts_ = new double [getNumCols()];
1527    CoinFillN(reducedCosts_,getNumCols(),0.);
1528  }
1529  return reducedCosts_;
1530}
1531
1532/** Get pointer to array[getNumRows()] of row activity levels (constraint
1533    matrix times the solution vector */
1534const double *
1535IpoptInterface::getRowActivity() const
1536{
1537  return problem_->g_sol();
1538}
1539
1540/** Get how many iterations it took to solve the problem (whatever
1541    "iteration" mean to the solver.
1542*/
1543int
1544IpoptInterface::getIterationCount() const
1545{
1546  if(IsValid(app_->Statistics()))
1547    return app_->Statistics()->IterationCount();
1548  else
1549    return 0;
1550}
1551
1552
1553/** Set a single column lower bound.
1554    Use -getInfinity() for -infinity. */
1555void
1556IpoptInterface::setColLower( int elementIndex, double elementValue )
1557{
1558  //  if(fabs(problem_->x_l()[elementIndex]-elementValue)>1e-06)
1559  problem_->SetVariableLowerBound(elementIndex,elementValue);
1560  hasBeenOptimized_ = false;
1561}
1562
1563/** Set a single column upper bound.
1564    Use getInfinity() for infinity. */
1565void
1566IpoptInterface::setColUpper( int elementIndex, double elementValue )
1567{
1568  //  if(fabs(problem_->x_u()[elementIndex]-elementValue)>1e-06)
1569  problem_->SetVariableUpperBound(elementIndex,elementValue);
1570  hasBeenOptimized_ = false;
1571}
1572
1573/** Set a single row lower bound.
1574    Use -getInfinity() for -infinity. */
1575void
1576IpoptInterface::setRowLower( int elementIndex, double elementValue )
1577{
1578  throw SimpleError("Not implemented yet but should be if necessary.",
1579      "setRowLower");
1580  hasBeenOptimized_ = false;
1581}
1582
1583/** Set a single row upper bound.
1584    Use getInfinity() for infinity. */
1585void
1586IpoptInterface::setRowUpper( int elementIndex, double elementValue )
1587{
1588  throw SimpleError("Not implemented yet but should be if necessary.",
1589      "setRowUpper");
1590  hasBeenOptimized_ = false;
1591}
1592
1593/** Set the type of a single row */
1594void
1595IpoptInterface::setRowType(int index, char sense, double rightHandSide,
1596    double range)
1597{
1598  throw SimpleError("Not implemented yet but should be if necessary.",
1599      "setRowType");
1600  hasBeenOptimized_ = false;
1601}
1602
1603
1604/// Set the objective function sense.
1605/// (1 for min (default), -1 for max)
1606void
1607IpoptInterface::setObjSense(double s)
1608{
1609  throw SimpleError("Can not change objective sense of an Ipopt problem.",
1610      "setObjSense");
1611  hasBeenOptimized_ = false;
1612}
1613
1614/** Set the primal solution variable values
1615
1616colsol[getNumCols()] is an array of values for the primal variables.
1617These values are copied to memory owned by the solver interface object
1618or the solver.  They will be returned as the result of getColSolution()
1619until changed by another call to setColSolution() or by a call to any
1620solver routine.  Whether the solver makes use of the solution in any
1621way is solver-dependent.
1622*/
1623void
1624IpoptInterface::setColSolution(const double *colsol)
1625{
1626  problem_->setxInit(getNumCols(), colsol);
1627  hasBeenOptimized_ = false;
1628  //  problem_->SetStartingPoint(getNumCols(), colsol);
1629}
1630
1631/** Set dual solution variable values
1632
1633rowprice[getNumRows()] is an array of values for the dual
1634variables. These values are copied to memory owned by the solver
1635interface object or the solver.  They will be returned as the result of
1636getRowPrice() until changed by another call to setRowPrice() or by a
1637call to any solver routine.  Whether the solver makes use of the
1638solution in any way is solver-dependent.
1639*/
1640
1641void
1642IpoptInterface::setRowPrice(const double * rowprice)
1643{
1644  problem_->setDualsInit(getNumCols()*2 + getNumRows(), rowprice);
1645  hasBeenOptimized_ = false;
1646}
1647
1648
1649/** Set the index-th variable to be a continuous variable */
1650void
1651IpoptInterface::setContinuous(int index)
1652{
1653  problem_->SetVariableType(index, Ipopt::TMINLP::CONTINUOUS);
1654  hasBeenOptimized_ = false;
1655}
1656/** Set the index-th variable to be an integer variable */
1657void
1658IpoptInterface::setInteger(int index)
1659{
1660  problem_->SetVariableType(index, Ipopt::TMINLP::INTEGER);
1661  hasBeenOptimized_ = false;
1662}
1663
1664/// Get objective function value (can't use default)
1665double
1666IpoptInterface::getObjValue() const
1667{
1668  return problem_->obj_value();
1669}
1670
1671
1672CoinWarmStart* IpoptInterface::getEmptyWarmStart () const
1673{
1674  return (dynamic_cast<CoinWarmStart *>(new IpoptWarmStart(1)));
1675}
1676
1677//#############################################################################
1678// Parameter related methods
1679//#############################################################################
1680
1681bool
1682IpoptInterface::setIntParam(OsiIntParam key, int value)
1683{
1684  //  debugMessage("OsiCpxSolverInterface::setIntParam(%d, %d)\n", key, value);
1685
1686  bool retval = false;
1687  switch (key) {
1688  case OsiMaxNumIteration:
1689    retval = false;
1690    break;
1691  case OsiMaxNumIterationHotStart:
1692    if( value >= 0 ) {
1693      retval = false;
1694    }
1695    else
1696      retval = false;
1697    break;
1698  case OsiLastIntParam:
1699    retval = false;
1700    break;
1701  }
1702  return retval;
1703}
1704
1705//-----------------------------------------------------------------------------
1706
1707bool
1708IpoptInterface::setDblParam(OsiDblParam key, double value)
1709{
1710  //  debugMessage("IpoptInterface::setDblParam(%d, %g)\n", key, value);
1711
1712  bool retval = false;
1713  switch (key) {
1714  case OsiDualObjectiveLimit:
1715    OsiDualObjectiveLimit_ = value;
1716    retval = true;
1717    break;
1718  case OsiPrimalObjectiveLimit:
1719    std::cerr<<"Can not set primal objective limit parameter"<<std::endl;
1720    retval = false;
1721    break;
1722  case OsiDualTolerance:
1723    std::cerr<<"Can not set dual tolerance parameter"<<std::endl;
1724    retval = false;
1725    break;
1726  case OsiPrimalTolerance:
1727    std::cerr<<"Can not set primal tolerance parameter"<<std::endl;
1728    retval = false;
1729  case OsiObjOffset:
1730    retval = OsiSolverInterface::setDblParam(key,value);
1731    break;
1732  case OsiLastDblParam:
1733    retval = false;
1734    break;
1735  }
1736  return retval;
1737}
1738
1739
1740//-----------------------------------------------------------------------------
1741
1742bool
1743IpoptInterface::setStrParam(OsiStrParam key, const std::string & value)
1744{
1745  //  debugMessage("IpoptInterface::setStrParam(%d, %s)\n", key, value.c_str());
1746
1747  bool retval=false;
1748  switch (key) {
1749  case OsiProbName:
1750    OsiSolverInterface::setStrParam(key,value);
1751    return retval = true;
1752  case OsiSolverName:
1753    return false;
1754  case OsiLastStrParam:
1755    return false;
1756  }
1757  return false;
1758}
1759
1760//-----------------------------------------------------------------------------
1761
1762bool
1763IpoptInterface::getIntParam(OsiIntParam key, int& value) const
1764{
1765  //  debugMessage("IpoptInterface::getIntParam(%d)\n", key);
1766
1767  bool retval = false;
1768  switch (key) {
1769  case OsiMaxNumIteration:
1770    retval = false;
1771    break;
1772  case OsiMaxNumIterationHotStart:
1773    retval = false;
1774    break;
1775  case OsiLastIntParam:
1776    retval = false;
1777    break;
1778  }
1779  return retval;
1780}
1781
1782//-----------------------------------------------------------------------------
1783
1784bool
1785IpoptInterface::getDblParam(OsiDblParam key, double& value) const
1786{
1787  //  debugMessage("IpoptInterface::getDblParam(%d)\n", key);
1788
1789  bool retval = false;
1790  switch (key) {
1791  case OsiDualObjectiveLimit:
1792    value = OsiDualObjectiveLimit_;
1793    retval = true;
1794    break;
1795  case OsiPrimalObjectiveLimit:
1796    value = getInfinity();
1797    retval = true;
1798    break;
1799  case OsiDualTolerance:
1800    retval = false;
1801    break;
1802  case OsiPrimalTolerance:
1803    retval = false;
1804    break;
1805  case OsiObjOffset:
1806    retval = OsiSolverInterface::getDblParam(key, value);
1807    break;
1808  case OsiLastDblParam:
1809    retval = false;
1810    break;
1811  }
1812  return retval;
1813}
1814
1815
1816//-----------------------------------------------------------------------------
1817
1818bool
1819IpoptInterface::getStrParam(OsiStrParam key, std::string & value) const
1820{
1821  //  debugMessage("IpoptInterface::getStrParam(%d)\n", key);
1822
1823  switch (key) {
1824  case OsiProbName:
1825    OsiSolverInterface::getStrParam(key, value);
1826    break;
1827  case OsiSolverName:
1828    value = "Ipopt";
1829    break;
1830  case OsiLastStrParam:
1831    return false;
1832  }
1833
1834  return true;
1835}
1836
1837void
1838IpoptInterface::turnOffIpoptOutput()
1839{
1840  std::string opt="print_level";
1841  app_->Options()->SetIntegerValue(opt, (Index)J_NONE,true);
1842}
1843void
1844IpoptInterface::turnOnIpoptOutput()
1845{
1846  std::string opt="print_level";
1847  app_->Options()->SetIntegerValue(opt, (Index)J_SUMMARY,true);
1848}
1849
1850void
1851IpoptInterface::randomStartingPoint()
1852{
1853  int numcols = getNumCols();
1854  const double * colLower = getColLower();
1855  const double * colUpper = getColUpper();
1856  double * sol = new double[numcols];
1857  for(int i = 0 ; i < numcols ; i++) {
1858    double lower = min(-maxRandomRadius_,colUpper[i] - maxRandomRadius_);
1859    lower = max(colLower[i], lower);
1860    double upper = max(maxRandomRadius_,colLower[i] + maxRandomRadius_);
1861    upper = min(colUpper[i],upper);
1862    lower = min(upper,lower);
1863    upper = max(upper, lower);
1864    double interval = upper - lower;
1865    sol[i] = CoinDrand48()*(interval) + lower;
1866//    printf("%f in [%f,%f]\n",sol[i],lower,upper);
1867    //  std::cout<<interval<<"\t";
1868  }
1869  //std::cout<<std::endl;
1870  unsetWarmStartOptions();
1871  setColSolution(sol);
1872  delete [] sol;
1873}
1874
1875/*******************************************************************************/
1876// Class for throwing errors reported from Ipopt
1877/******************************************************************************/
1878
1879std::string
1880IpoptInterface::UnsolvedError::errorNames[17] ={"Solve succeeded",
1881    "Solved to acceptable level",
1882    "Infeasible Problem Detected",
1883    "Search_Direction_Becomes_Too_Small",
1884    "Diverging iterates",
1885    "User Requested Stop",
1886    "Maximum Iterations Exceeded",
1887    "Restoration Failed",
1888    "Error in step computation",
1889    "Not enough degrees of freedom",
1890    "Invalid problem definition",
1891    "Invalid option",
1892    "Invalid number detected",
1893    "Unrecoverable exception",
1894    "NonIpopt exception thrown",
1895    "Insufficient memory",
1896    "Internal error"};
1897
1898IpoptInterface::UnsolvedError::UnsolvedError(int errorNum)
1899{
1900  errorNum_ = errorNum;
1901}
1902
1903const std::string &
1904IpoptInterface::UnsolvedError::errorName() const
1905{
1906  if(errorNum_ >=0)
1907    return errorNames[errorNum_];
1908  if(errorNum_ == -1) return errorNames[6];
1909  else if(errorNum_ == -2) return errorNames[7];
1910  else if(errorNum_ == -3) return errorNames[8];
1911  else if(errorNum_ == -10) return errorNames[9];
1912  else if(errorNum_ == -11) return errorNames[10];
1913  else if(errorNum_ == -12) return errorNames[11];
1914  else if(errorNum_ == -13) return errorNames[12];
1915  else if(errorNum_ == -100) return errorNames[13];
1916  else if(errorNum_ == -101) return errorNames[14];
1917  else if(errorNum_ == -102) return errorNames[15];
1918  else if(errorNum_ == -199) return errorNames[16];
1919  throw IpoptInterface::SimpleError("UnsolvedError::errorName()","Unrecognized optimization status in ipopt.");
1920}
1921
1922void
1923IpoptInterface::UnsolvedError::printError(std::ostream &os)
1924{
1925  os<<"Ipopt exited with error code "<<errorNum_<<" "<<errorName()<<std::endl;
1926}
1927
1928
1929/** This methods initialiaze arrays for storing the jacobian */
1930int IpoptInterface::initializeJacobianArrays()
1931{
1932  Ipopt::Index n, m, nnz_h_lag;
1933  Ipopt::TNLP::IndexStyleEnum index_style;
1934  tminlp_->get_nlp_info( n, m, nnz_jac, nnz_h_lag, index_style);
1935
1936  if(jRow_ != NULL) delete jRow_;
1937  if(jCol_ != NULL) delete jCol_;
1938  if(jValues_ != NULL) delete jValues_;
1939
1940  jRow_ = new Ipopt::Index[nnz_jac];
1941  jCol_ = new Ipopt::Index[nnz_jac];
1942  jValues_ = new Ipopt::Number[nnz_jac];
1943  tminlp_->eval_jac_g(n, NULL, 0, m, nnz_jac, jRow_, jCol_, NULL);
1944
1945
1946  if(constTypes_ != NULL) delete [] constTypes_;
1947//  if(constTypesNum_ != NULL) delete [] constTypesNum_;
1948
1949  constTypes_ = new Ipopt::TMINLP::ConstraintType[getNumRows()];
1950  tminlp_->get_constraints_types(getNumRows(), constTypes_);
1951//  constTypesNum_ = new int[getNumRows()];
1952  for(int i = 0; i < getNumRows() ; i++) {
1953    if(constTypes_[i]==Ipopt::TMINLP::LINEAR) {
1954      //constTypesNum_[i] =
1955      nLinear_++;
1956    }
1957    else if(constTypes_[i]==Ipopt::TMINLP::NON_LINEAR) {
1958      //constTypesNum_[i] =
1959      nNonLinear_++;
1960    }
1961  }
1962  return nnz_jac;
1963}
1964
1965
1966/** get NLP constraint violation of current point */
1967double
1968IpoptInterface::getConstraintViolation()
1969{
1970  int numcols = getNumCols();
1971  int numrows = getNumRows();
1972  double * g = new double[numrows];
1973  tminlp_->eval_g(numcols,getColSolution(),1,numrows,g);
1974  const double * rowLower = getRowLower();
1975  const double * rowUpper = getRowUpper();
1976
1977
1978  double norm = 0;
1979  for(int i = 0; i< numrows ; i++) {
1980    if(constTypes_[i] == Ipopt::TMINLP::NON_LINEAR) {
1981      double rowViolation = max(0.,rowLower[i] - g[i]);
1982      if(rowViolation  > 0.) rowViolation  /= fabs(rowLower[i]);
1983      double rowViolation2 = max(0.,g[i] - rowUpper[i]);
1984      if(rowViolation2 > 0.) rowViolation2 /= fabs(rowUpper[i]);
1985      norm = max(rowViolation, norm);
1986    }
1987  }
1988  return norm;
1989}
1990
1991
1992//A procedure to try to remove small coefficients in OA cuts (or make it non small
1993static inline
1994bool cleanNnz(double &value, double colLower, double colUpper,
1995    double rowLower, double rowUpper, double colsol,
1996    double & lb, double &ub, double tiny, double veryTiny)
1997{
1998  if(fabs(value)>= tiny) return 1;
1999
2000  if(fabs(value)<veryTiny) return 0;//Take the risk?
2001
2002  //try and remove
2003  double infty = 1e20;
2004  bool colUpBounded = colUpper < 10000;
2005  bool colLoBounded = colLower > -10000;
2006  bool rowNotLoBounded =  rowLower <= - infty;
2007  bool rowNotUpBounded = rowUpper >= infty;
2008  bool pos =  value > 0;
2009
2010  if(!rowNotLoBounded && ! rowNotUpBounded)//would have to either choose side or duplicate cut
2011  {
2012    std::cerr<<"OA on non-convex constraint is very experimental, don't know how to remove"
2013    <<std::endl;
2014  }
2015
2016  if(colLoBounded && pos && rowNotUpBounded) {
2017    lb += value * (colsol - colLower);
2018    return 0;
2019  }
2020  else
2021    if(colLoBounded && !pos && rowNotLoBounded) {
2022      ub += value * (colsol - colLower);
2023      return 0;
2024    }
2025    else
2026      if(colUpBounded && !pos && rowNotUpBounded) {
2027        lb += value * (colsol - colUpper);
2028        return 0;
2029      }
2030      else
2031        if(colUpBounded && pos && rowNotLoBounded) {
2032          ub += value * (colsol - colUpper);
2033          return 0;
2034        }
2035  //can not remove coefficient increase it to smallest non zero
2036  if(pos) value = tiny;
2037  else
2038    value = - tiny;
2039  return 1;
2040}
2041
2042/** Get the outer approximation constraints at the current optimum of the
2043  current ipopt problem */
2044void
2045IpoptInterface::getOuterApproximation(OsiCuts &cs, bool getObj)
2046{
2047  int n,m, nnz_jac_g, nnz_h_lag;
2048  Ipopt::TNLP::IndexStyleEnum index_style;
2049  tminlp_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2050  if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2051    initializeJacobianArrays();
2052  assert(jRow_ != NULL);
2053  assert(jCol_ != NULL);
2054  double * g = new double[m];
2055  tminlp_->eval_jac_g(n, getColSolution(), 1, m, nnz_jac_g, NULL, NULL, jValues_);
2056  tminlp_->eval_g(n,getColSolution(),1,m,g);
2057  //As jacobian is stored by cols fill OsiCuts with cuts
2058  CoinPackedVector * cuts = new CoinPackedVector[nNonLinear_ + 1];
2059  double * lb = new double[nNonLinear_ + 1];
2060  double * ub = new double[nNonLinear_ + 1];
2061  int * binding = new int[nNonLinear_ + 1];//store binding constraints at current opt (which are added to OA) -1 if constraint is not binding, otherwise index in subset of binding constraints
2062  int numBindings = 0;
2063   
2064  const double * rowLower = getRowLower();
2065  const double * rowUpper = getRowUpper();
2066  const double * colLower = getColLower();
2067  const double * colUpper = getColUpper();
2068  const double * duals = getRowPrice();
2069  double infty = 1e100 *getInfinity();
2070 
2071  for(int i = 0; i< m ; i++) {
2072    if(constTypes_[i] == Ipopt::TMINLP::NON_LINEAR) {
2073      if(rowLower[i] > -infty && rowUpper[i] < infty && fabs(duals[i]) == 0.)
2074      {
2075        binding[i] = -1;
2076        std::cerr<<"non binding constraint"<<std::endl;
2077        continue;
2078      }
2079      binding[i] = numBindings;
2080      if(rowLower[i] > - infty_)
2081        lb[numBindings] = rowLower[i] - g[i];
2082      else
2083        lb[numBindings] = - infty;
2084      if(rowUpper[i] < infty_)
2085        ub[numBindings] = rowUpper[i] - g[i];
2086      else
2087        ub[numBindings] = infty_;
2088      if(rowLower[i] > -infty && rowUpper[i] < infty)
2089      {
2090        if(duals[i] >= 0)// <= inequality
2091          lb[numBindings] = - infty;
2092        if(duals[i] <= 0)// >= inequality
2093          ub[numBindings] = infty;
2094      }
2095     
2096      numBindings++;
2097    }
2098  }
2099
2100  for(int i = 0 ; i < nnz_jac_g ; i++) {
2101    if(constTypes_[jRow_[i] - 1] == Ipopt::TMINLP::NON_LINEAR) {
2102      //"clean" coefficient
2103      if(cleanNnz(jValues_[i],colLower[jCol_[i] - 1], colUpper[jCol_[i]-1],
2104          rowLower[jRow_[i] - 1], rowUpper[jRow_[i] - 1],
2105          getColSolution()[jCol_[i] - 1],
2106          lb[binding[jRow_[i] - 1]],
2107          ub[binding[jRow_[i] - 1]], tiny_, veryTiny_)) {
2108        cuts[binding[jRow_[i] - 1]].insert(jCol_[i]-1,jValues_[i]);
2109        if(lb[binding[jRow_[i] - 1]] > - infty_)
2110          lb[binding[jRow_[i] - 1]] += jValues_[i] * getColSolution()[jCol_ [i] - 1];
2111        if(ub[binding[jRow_[i] - 1]] < infty_)
2112        ub[binding[jRow_[i] - 1]] += jValues_[i] * getColSolution()[jCol_ [i] - 1];
2113      }
2114    }
2115  }
2116
2117  for(int i = 0; i< numBindings ; i++) {
2118    OsiRowCut newCut;
2119    //    if(lb[i]>-1e20) assert (ub[i]>1e20);
2120
2121    newCut.setGloballyValid();
2122    newCut.setLb(lb[i]);
2123    newCut.setUb(ub[i]);
2124    newCut.setRow(cuts[i]);
2125    //    CftValidator validator;
2126    //    validator(newCut);
2127    cs.insert(newCut);
2128  }
2129
2130  delete[] g;
2131  delete [] cuts;
2132
2133  if(getObj)  // Get the objective cuts
2134  {
2135    double * obj = new double [n];
2136    tminlp_->eval_grad_f(n, getColSolution(), 1,obj);
2137    double f;
2138    tminlp_->eval_f(n, getColSolution(), 1, f);
2139
2140    CoinPackedVector v;
2141    v.reserve(n);
2142    lb[nNonLinear_] = -f;
2143    ub[nNonLinear_] = -f;
2144    //double minCoeff = 1e50;
2145    for(int i = 0; i<n ; i++)
2146    {
2147      if(cleanNnz(obj[i],colLower[i], colUpper[i],
2148          -getInfinity(), 0,
2149          getColSolution()[i],
2150          lb[nNonLinear_],
2151          ub[nNonLinear_],tiny_, 1e-15)) {
2152        //            minCoeff = min(fabs(obj[i]), minCoeff);
2153        v.insert(i,obj[i]);
2154        lb[nNonLinear_] += obj[i] * getColSolution()[i];
2155        ub[nNonLinear_] += obj[i] * getColSolution()[i];
2156      }
2157    }
2158    v.insert(n,-1);
2159    OsiRowCut newCut;
2160    newCut.setGloballyValid();
2161    newCut.setRow(v);
2162    newCut.setLb(-DBL_MAX/*Infinity*/);
2163    newCut.setUb(ub[nNonLinear_]);
2164//     CftValidator validator;
2165//     validator(newCut);
2166    cs.insert(newCut);
2167    delete [] obj;
2168  }
2169
2170  delete []lb;
2171  delete[]ub;
2172}
2173
2174double
2175IpoptInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs)
2176{
2177  if(! IsValid(feasibilityProblem_)) {
2178    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2179  }
2180  feasibilityProblem_->set_dist2point_obj(n,(const Ipopt::Number *) x_bar,(const Ipopt::Index *) inds);
2181  nCallOptimizeTNLP_++;
2182  totalNlpSolveTime_-=CoinCpuTime();
2183  Ipopt::IpoptApplication app2;
2184  app2.Options()->SetIntegerValue("print_level", (Ipopt::Index) 0);
2185  optimization_status_ = app2.OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2186  totalNlpSolveTime_+=CoinCpuTime();
2187  getOuterApproximation(cs, 0);
2188  hasBeenOptimized_=true;
2189  return getObjValue();
2190}
2191
2192
2193
2194
2195
2196void
2197IpoptInterface::extractLinearRelaxation(OsiSolverInterface &si, bool getObj)
2198{
2199  initialSolve();
2200
2201  CoinPackedMatrix mat;
2202  double * rowLow = NULL;
2203  double * rowUp = NULL;
2204
2205  int n;
2206  int m;
2207  int nnz_jac_g;
2208  int nnz_h_lag;
2209  Ipopt::TNLP::IndexStyleEnum index_style;
2210  //Get problem information
2211  tminlp_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2212
2213  //if not allocated allocate spaced for stroring jacobian
2214  if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2215    initializeJacobianArrays();
2216
2217  //get Jacobian
2218  tminlp_->eval_jac_g(n, getColSolution(), 1, m, nnz_jac_g, NULL, NULL, jValues_);
2219
2220
2221  double *g = new double[m];
2222  tminlp_->eval_g(n, getColSolution(),1,m,g);
2223
2224  rowLow = new double[m];
2225  rowUp = new double[m];
2226  int * binding = new int[m];//store binding constraints at current opt (which are added to OA) -1 if constraint is not binding, otherwise index in subset of binding constraints
2227  int numBindings = 0;
2228  const double * rowLower = getRowLower();
2229  const double * rowUpper = getRowUpper();
2230  const double * colLower = getColLower();
2231  const double * colUpper = getColUpper();
2232  const double * duals = getRowPrice();
2233  assert(m==getNumRows());
2234  double infty = getInfinity() * 1e100;
2235  for(int i = 0 ; i < m ; i++) {
2236    {
2237      if(constTypes_[i] == Ipopt::TMINLP::NON_LINEAR) {
2238        //If constraint is equality not binding do not add
2239        if(rowLower[i] > -infty && rowUpper[i] < infty && fabs(duals[i]) == 0.)
2240        {
2241            binding[i] = -1;
2242            continue;
2243        }
2244        else
2245          binding[i] = numBindings;
2246        if(rowLower[i] > - infty_)
2247          rowLow[numBindings] = (rowLower[i] - g[i]) - 1e-07;
2248        else
2249          rowLow[numBindings] = - infty_;
2250        if(rowUpper[i] < infty_)
2251          rowUp[numBindings] =  (rowUpper[i] - g[i]) + 1e-07;
2252        else
2253          rowUp[numBindings] = infty_;
2254       
2255        //If equality or ranged constraint only add one side by looking at sign of dual multiplier
2256        if(rowLower[i] > -infty && rowUpper[i] < infty)
2257        {
2258          if(duals[i] >= 0)// <= inequality
2259            rowLow[numBindings] = - infty;
2260          if(duals[i] <= 0)// >= inequality
2261            rowUp[numBindings] = infty;
2262        }
2263        numBindings++;
2264      }
2265      else {
2266        binding[i] = numBindings;//All are considered as bindings
2267        rowLow[numBindings] = (rowLower[i] - g[i]);
2268        rowUp[numBindings] =  (rowUpper[i] - g[i]);
2269        numBindings++;
2270      }
2271    }
2272  }
2273
2274  //Then convert everything to a CoinPackedMatrix (col ordered)
2275  CoinBigIndex * inds = new CoinBigIndex[nnz_jac_g + 1];
2276  double * vals = new double [nnz_jac_g + 1];
2277  int * start = new int[n+1];
2278  int * length = new int[n];
2279  bool needOrder = false;
2280  int nnz = 0;
2281  if(nnz_jac_g > 0)
2282  {
2283  for(int k = 0; k < jCol_[0]; k++) {
2284    start[k] = nnz;
2285    length[k] = 0;
2286  }
2287  for(int i = 0 ; i < nnz_jac_g - 1 ; i++) {
2288    {
2289      if(jCol_[i + 1] < jCol_ [i] || ( jCol_ [i+1] == jCol_[i] && jRow_[i + 1] <= jRow_[i]) ) {
2290        needOrder = 1;
2291        break;
2292      }
2293      if(Ipopt::TMINLP::LINEAR //Always accept coefficients from linear constraints
2294          || //For other clean tinys
2295          cleanNnz(jValues_[i],colLower[jCol_[i] - 1], colUpper[jCol_[i]-1],
2296              rowLower[jRow_[i] - 1], rowUpper[jRow_[i] - 1],
2297              getColSolution()[jCol_[i] - 1],
2298              rowLow[binding[jRow_[i] - 1]],
2299              rowUp[binding[jRow_[i] -1]], tiny_, veryTiny_)) {
2300        if(binding[jRow_[i] - 1] == -1) continue;
2301        vals[nnz] = jValues_[i];
2302        if(rowLow[binding[jRow_[i] - 1]] > - infty_)
2303          rowLow[binding[jRow_[i] - 1]] += jValues_[i] * getColSolution()[jCol_ [i] - 1];
2304        if(rowUp[binding[jRow_[i] - 1]] < infty_)
2305          rowUp[binding[jRow_[i] -1]] += jValues_[i] *getColSolution()[jCol_[i] -1];
2306        inds[nnz] = binding[jRow_[i] - 1];//jRow_[i ] - 1;
2307        length[jCol_[i] - 1]++;
2308        nnz++;
2309      }
2310    }
2311    if(jCol_[i + 1] > jCol_[i]) {
2312      for(int k = jCol_[i]; k < jCol_[i + 1]; k++) {
2313        start[k] = nnz;
2314        length[k] = 0;
2315      }
2316    }
2317  }
2318  if(!needOrder) {
2319    {
2320      length[jCol_[nnz_jac_g -1] - 1]++;
2321      vals[nnz] = jValues_[nnz_jac_g - 1];
2322      rowLow[binding[jRow_[nnz_jac_g - 1] - 1]] += jValues_[nnz_jac_g - 1] * getColSolution()[jCol_ [nnz_jac_g - 1] - 1];
2323      rowUp[binding[jRow_[nnz_jac_g - 1] -1]] += jValues_[nnz_jac_g - 1] *getColSolution()[jCol_[nnz_jac_g - 1] -1];
2324      inds[nnz++] = binding[jRow_[nnz_jac_g - 1] - 1];
2325    }
2326    for(int i = jCol_[nnz_jac_g -1] ; i < n ; i++) {
2327      start[i] = nnz;
2328      length[i] = 0;
2329    }
2330    start[n]=nnz;
2331  }
2332  else {
2333    std::cerr<<"jacobian matrix is not ordered properly"<<std::endl;
2334    throw -1;
2335  }
2336  delete [] g;
2337  }
2338  else {
2339   for (int i = 0 ; i < n ; i++)
2340   {
2341     length[i] = 0;
2342     start[i] = 0;
2343   }
2344   start[n]=0;
2345 }
2346 mat.assignMatrix(false, m, n, nnz, vals, inds, start, length);
2347  int numcols=getNumCols();
2348  double *obj = new double[numcols];
2349  for(int i = 0 ; i < numcols ; i++)
2350    obj[i] = 0;
2351  mat.transpose();
2352 
2353#if 0
2354  std::cout<<"Mat is ordered by "<<
2355  <<mat.isColOrdered?"columns.":"rows."
2356  <<std::endl;
2357#endif
2358 
2359  si.loadProblem(mat, getColLower(), getColUpper(), obj, rowLow, rowUp);
2360  delete [] rowLow;
2361  delete [] rowUp;
2362  delete [] binding;
2363  for(int i = 0 ; i < getNumCols() ; i++) {
2364    if(isInteger(i))
2365      si.setInteger(i);
2366  }
2367  if(getObj) {
2368    //add variable alpha
2369    //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
2370    CoinPackedVector a;
2371    si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
2372
2373    // Now get the objective cuts
2374    // get the gradient, pack it and add the cut
2375    tminlp_->eval_grad_f(n, getColSolution(), 1,obj);
2376    double ub;
2377    tminlp_->eval_f(n, getColSolution(), 1, ub);
2378    ub*=-1;
2379    double lb = -1e300;
2380    CoinPackedVector objCut;
2381    CoinPackedVector * v = &objCut;
2382    v->reserve(n);
2383    for(int i = 0; i<n ; i++) {
2384     if(nnz_jac_g)
2385     {
2386      if(cleanNnz(obj[i],colLower[i], colUpper[i],
2387          -getInfinity(), 0,
2388          getColSolution()[i],
2389          lb,
2390          ub, tiny_, veryTiny_)) {
2391        v->insert(i,obj[i]);
2392        lb += obj[i] * getColSolution()[i];
2393        ub += obj[i] * getColSolution()[i];
2394      }
2395     }
2396     else //Unconstrained problem can not put clean coefficient
2397     {
2398         if(cleanNnz(obj[i],colLower[i], colUpper[i],
2399          -getInfinity(), 0,
2400          getColSolution()[i],
2401          lb,
2402          ub, 1e-03, 1e-08)) {
2403        v->insert(i,obj[i]);
2404        lb += obj[i] * getColSolution()[i];
2405        ub += obj[i] * getColSolution()[i];
2406         }
2407     }
2408    }
2409    v->insert(n,-1);
2410    si.addRow(objCut, lb, ub);
2411    delete [] obj;
2412  }
2413
2414si.writeMps("init");
2415
2416  setWarmStartOptions();
2417  setColSolution(problem()->x_sol());
2418  setRowPrice(problem()->duals_sol());
2419
2420}
2421
2422
Note: See TracBrowser for help on using the repository browser.