source: trunk/Bonmin/src/IpoptInterface/IpoptInterface.cpp @ 1

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

imported initial code

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 70.4 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,5,1,
45      "Set the level of output of the branch-and-bound : "
46      "0 - none, 1 - minimal, 2 - normal low, 3 - normal high, 4 - verbose"
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,5,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,5,0,
65      "Set the level of output of the MILP subsolver in OA : "
66      "0 - none, 1 - minimal, 2 - normal low, 3 - normal high, 4 - verbose"
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{}
433
434
435/** Constructor with given IpSolver and TMINLP */
436IpoptInterface::IpoptInterface (Ipopt::SmartPtr<Ipopt::TMINLP> tminlp):
437    OsiSolverInterface(),
438    tminlp_(tminlp),
439    problem_(NULL),
440    app_(NULL),
441    rowsense_(NULL),
442    rhs_(NULL),
443    rowrange_(NULL),
444    reducedCosts_(NULL),
445    OsiDualObjectiveLimit_(1e200),
446//    optimization_status_(HasNotBeenOptimized),
447    varNames_(NULL),
448    hasVarNamesFile_(true),
449    nCallOptimizeTNLP_(0),
450    totalNlpSolveTime_(0),
451    totalIterations_(0),
452    maxRandomRadius_(1e08),
453    pushValue_(1e-02),
454    numRetryInitial_(-1),
455    numRetryResolve_(-1),
456    numRetryUnsolved_(1),
457    ipoptIMessages_(),
458    pretendFailIsInfeasible_(false),
459    hasContinuedAfterNlpFailure_(false),
460    numIterationSuspect_(-1),
461    hasBeenOptimized_(false),
462    warmStartStrategy_(1),
463    obj_(NULL),
464    feasibilityProblem_(NULL),
465    jRow_(NULL),
466    jCol_(NULL),
467    jValues_(NULL),
468    nnz_jac(0),
469    constTypes_(NULL),
470    constTypesNum_(NULL),
471    nLinear_(0),
472    nNonLinear_(0),
473    tiny_(1e-08),
474    veryTiny_(1e-17)
475{
476  assert(IsValid(tminlp));
477  app_ = new Ipopt::IpoptApplication();
478
479  SmartPtr<RegisteredOptions> roptions = app_->RegOptions();
480  register_ALL_options(roptions);
481  app_->Initialize("");
482  extractInterfaceParams();
483  set_ipopt_minlp_default(app_->Options());
484  // set the default options... expect_infeasible, etc...
485  problem_ = new Ipopt::TMINLP2TNLP(tminlp_, *app_->Options());
486  feasibilityProblem_ =
487    new Ipopt::TNLP2FPNLP
488    (Ipopt::SmartPtr<Ipopt::TNLP>(Ipopt::GetRawPtr(problem_)));
489
490}
491
492void
493IpoptInterface::readOptionFile(const char * fileName)
494{
495  app_->Initialize(fileName);
496  extractInterfaceParams();
497  set_ipopt_minlp_default(app_->Options());
498}
499void
500IpoptInterface::extractInterfaceParams()
501{
502  if (IsValid(app_)) {
503    app_->Options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,"bonmin.");
504    app_->Options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,"bonmin.");
505    app_->Options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,"bonmin.");
506    app_->Options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,"bonmin.");
507    app_->Options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,"bonmin.");
508    app_->Options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,"bonmin.");
509    app_->Options()->GetEnumValue("warm_start",warmStartStrategy_,"bonmin.");
510    app_->Options()->GetNumericValue
511    ("warm_start_bound_frac" ,pushValue_,"bonmin.");
512    app_->Options()->GetNumericValue("tiny_element",tiny_,"bonmin.");
513    app_->Options()->GetNumericValue("very_tiny_element",veryTiny_,"bonmin.");
514  }
515}
516
517/// Clone
518OsiSolverInterface * IpoptInterface::clone(bool CopyData) const
519{
520  //    debugMessage("IpoptInterface::clone(%d)\n", CopyData);
521  static int debug_no_clone = 0;
522  debug_no_clone++;
523  //    assert(debug_no_clone==1);
524  return new IpoptInterface(*this);
525}
526
527/// Copy constructor
528IpoptInterface::IpoptInterface (const IpoptInterface &source):
529    OsiSolverInterface(source),
530    tminlp_(source.tminlp_),
531    problem_(NULL),
532    app_(NULL),
533    rowsense_(NULL),
534    rhs_(NULL),
535    rowrange_(NULL),
536    reducedCosts_(NULL),
537    OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
538    optimization_status_(source.optimization_status_),
539    varNames_(NULL),
540    hasVarNamesFile_(source.hasVarNamesFile_),
541    nCallOptimizeTNLP_(0),
542    totalNlpSolveTime_(0),
543    totalIterations_(0),
544    maxRandomRadius_(source.maxRandomRadius_),
545    pushValue_(source.pushValue_),
546    numRetryInitial_(source.numRetryInitial_),
547    numRetryResolve_(source.numRetryResolve_),
548    numRetryUnsolved_(source.numRetryUnsolved_),
549    ipoptIMessages_(),
550    pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
551    hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
552    numIterationSuspect_(source.numIterationSuspect_),
553    hasBeenOptimized_(source.hasBeenOptimized_),
554    warmStartStrategy_(source.warmStartStrategy_),
555    obj_(NULL),
556    feasibilityProblem_(NULL),
557    jRow_(NULL),
558    jCol_(NULL),
559    jValues_(NULL),
560    nnz_jac(source.nnz_jac),
561    constTypes_(NULL),
562    constTypesNum_(NULL),
563    nLinear_(0),
564    nNonLinear_(0),
565    tiny_(source.tiny_),
566    veryTiny_(source.veryTiny_)
567{
568  // Copy options from old application
569  if(IsValid(source.tminlp_)) {
570    app_ = new Ipopt::IpoptApplication();
571
572    SmartPtr<RegisteredOptions> roptions = app_->RegOptions();
573    register_ALL_options(roptions);
574    // Copy the options
575    *app_->Options()=*source.app_->Options();
576
577    //    extractInterfaceParams();
578
579    problem_ = new Ipopt::TMINLP2TNLP(tminlp_, *app_->Options());
580    problem_->copyUserModification(*source.problem_);
581    pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
582  }
583  else {
584    throw SimpleError("Don't know how to copy an empty IpoptInterface.",
585        "copy constructor");
586  }
587
588  if(source.obj_) {
589    obj_ = new double[source.getNumCols()];
590    CoinCopyN(source.obj_, source.getNumCols(), obj_);
591  }
592  if(IsValid(source.tminlp_))
593    feasibilityProblem_ = new Ipopt::TNLP2FPNLP
594        (Ipopt::SmartPtr<Ipopt::TNLP>(Ipopt::GetRawPtr(problem_)));
595  else
596    throw CoinError("Don't know how to copy an empty IpoptOAInterface.",
597        "copy constructor","IpoptOAInterface");
598
599
600  if(source.jValues_!=NULL && source.jRow_ != NULL && source.jCol_ != NULL && nnz_jac>0) {
601    jValues_ = new double [nnz_jac];
602    jCol_    = new Ipopt::Index [nnz_jac];
603    jRow_    = new Ipopt::Index [nnz_jac];
604    CoinCopyN(source.jValues_ , nnz_jac,jValues_ );
605    CoinCopyN(source.jCol_    , nnz_jac,jCol_    );
606    CoinCopyN(source.jRow_    , nnz_jac,jRow_    );
607
608    if(source.constTypes_ != NULL) {
609      constTypes_ = new Ipopt::TMINLP::ConstraintType[getNumRows()];
610      CoinCopyN(source.constTypes_, getNumRows(), constTypes_);
611    }
612    if(source.constTypesNum_ != NULL) {
613      constTypesNum_ = new int[getNumRows()];
614      CoinCopyN(source.constTypesNum_, getNumRows(), constTypesNum_);
615    }
616  }
617  else if(nnz_jac > 0) {
618    throw CoinError("Arrays for storing jacobian are inconsistant.",
619        "copy constructor","IpoptOAInterface");
620  }
621  // Process the output options
622  app_->Initialize("");
623}
624
625
626Ipopt::SmartPtr<Ipopt::OptionsList> IpoptInterface::retrieve_options()
627{
628  if(!IsValid(app_)) {
629    std::cout<<"Can not parse options when no IpApplication has been created"<<std::endl;
630    return NULL;
631  }
632  else
633    return app_->Options();
634}
635
636/// Assignment operator
637IpoptInterface & IpoptInterface::operator=(const IpoptInterface& rhs)
638{
639  if(this!= &rhs) {
640    OsiSolverInterface::operator=(rhs);
641    OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
642    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
643    totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
644    totalIterations_ = rhs.totalIterations_;
645    maxRandomRadius_ = rhs.maxRandomRadius_;
646    hasVarNamesFile_ = rhs.hasVarNamesFile_;
647    pushValue_ = rhs.pushValue_;
648    optimization_status_ = rhs.optimization_status_;
649
650    if(IsValid(rhs.tminlp_)) {
651
652      tminlp_ = rhs.tminlp_;
653      app_ = new Ipopt::IpoptApplication();
654      problem_ = new Ipopt::TMINLP2TNLP(tminlp_, *app_->Options());
655
656      feasibilityProblem_ = new Ipopt::TNLP2FPNLP
657          (Ipopt::SmartPtr<Ipopt::TNLP>(Ipopt::GetRawPtr(problem_)));
658      nnz_jac = rhs.nnz_jac;
659
660      if(constTypes_ != NULL) {
661        delete [] constTypes_;
662        constTypes_ = NULL;
663      }
664      if(rhs.constTypes_ != NULL) {
665        constTypes_ = new Ipopt::TMINLP::ConstraintType[getNumRows()];
666        CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
667      }
668
669      if(constTypesNum_ != NULL) {
670        delete [] constTypesNum_;
671        constTypesNum_ = NULL;
672      }
673      if(rhs.constTypesNum_ != NULL) {
674        constTypesNum_ = new int[getNumRows()];
675        CoinCopyN(rhs.constTypesNum_, getNumRows(), constTypesNum_);
676      }
677
678      if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
679        jValues_ = new double [nnz_jac];
680        jCol_    = new Ipopt::Index [nnz_jac];
681        jRow_    = new Ipopt::Index [nnz_jac];
682        CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
683        CoinCopyN(rhs.jCol_    , nnz_jac,jCol_    );
684        CoinCopyN(rhs.jRow_    , nnz_jac,jRow_    );
685      }
686      else if(nnz_jac > 0) {
687        throw CoinError("Arrays for storing jacobian are inconsistant.",
688            "copy constructor",
689            "IpoptOAInterface");
690      }
691      tiny_ = rhs.tiny_;
692      veryTiny_ = rhs.veryTiny_;
693
694
695    }
696    else {
697      tminlp_ =NULL;
698      app_ = NULL;
699      problem_ = NULL;
700      feasibilityProblem_ = NULL;
701    }
702
703
704    if(obj_) {
705      delete [] obj_;
706      obj_ = NULL;
707    }
708    if(rhs.obj_) {
709      obj_ = new double[rhs.getNumCols()];
710      CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
711    }
712
713    delete [] varNames_;
714    varNames_ = NULL;
715
716    if(rhs.varNames_) {
717      rhs.varNames_ = new std::string[getNumCols()];
718      CoinCopyN(rhs.varNames_, getNumCols(), varNames_);
719    }
720
721    hasVarNamesFile_ = rhs.hasVarNamesFile_;
722
723    nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
724    totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
725    totalIterations_ = rhs.totalIterations_;
726    maxRandomRadius_ = rhs.maxRandomRadius_;
727    pushValue_ = rhs.pushValue_;
728    numRetryInitial_ = rhs.numRetryInitial_;
729    numRetryResolve_ = rhs.numRetryResolve_;
730    numRetryUnsolved_ = rhs.numRetryUnsolved_;
731    pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
732    numIterationSuspect_ = rhs.numIterationSuspect_;
733
734    hasBeenOptimized_ = rhs.hasBeenOptimized_;
735
736    freeCachedData();
737  }
738  return *this;
739}
740
741/// Destructor
742IpoptInterface::~IpoptInterface ()
743{
744  freeCachedData();
745  delete [] jRow_;
746  delete [] jCol_;
747  delete [] jValues_;
748  delete [] constTypes_;
749  delete [] constTypesNum_;
750}
751
752void
753IpoptInterface::freeCachedColRim()
754{
755  if(varNames_!=NULL) {
756    delete [] varNames_;
757    varNames_ = NULL;
758  }
759  if(reducedCosts_!=NULL) {
760    delete []  reducedCosts_;
761    reducedCosts_ = NULL;
762  }
763
764}
765
766void
767IpoptInterface::freeCachedRowRim()
768{
769  if(rowsense_!=NULL) {
770    delete []  rowsense_;
771    rowsense_ = NULL;
772  }
773  if(rhs_!=NULL) {
774    delete []  rhs_;
775    rhs_ = NULL;
776  }
777  if(rowrange_!=NULL) {
778    delete []  rowrange_;
779    rowrange_ = NULL;
780  }
781  //   if(dualsol_!=NULL)
782  //     {
783  //       delete []  dualsol_; dualsol_ = NULL;
784  //     }
785}
786
787void
788IpoptInterface::freeCachedData()
789{
790  freeCachedColRim();
791  freeCachedRowRim();
792}
793
794///////////////////////////////////////////////////////////////////
795// WarmStart Information                                                                           //
796///////////////////////////////////////////////////////////////////
797
798/// Get warmstarting information
799CoinWarmStart*
800IpoptInterface::getWarmStart() const
801{
802  if(warmStartStrategy_) {
803    if(warmStartStrategy_==2) {
804      SmartPtr<IpoptInteriorWarmStarter> warm_starter =
805        SmartPtr<IpoptInteriorWarmStarter>(problem_->GetWarmStarter());
806      return new IpoptWarmStart(*this, warm_starter);
807    }
808    else  return new IpoptWarmStart(*this, NULL);
809  }
810  else
811    return new IpoptWarmStart(getNumCols(), getNumRows());
812}
813
814
815bool
816IpoptInterface::setWarmStart(const CoinWarmStart* warmstart)
817{
818  if(!warmstart || !warmStartStrategy_)
819    return 0;
820  hasBeenOptimized_ = false;
821  const IpoptWarmStart * ws = dynamic_cast<const IpoptWarmStart*> (warmstart);
822  if(ws->empty())//reset initial point and leave
823  {
824    unsetWarmStartOptions();
825    return 1;
826  }
827  setWarmStartOptions();
828  int numcols = getNumCols();
829  int numrows = getNumRows();
830  const double * colLow = getColLower();
831  const double * colUp = getColUpper();
832  for(int i = 0; i < numcols ; i++) {
833    CoinWarmStartBasis::Status status = ws->getStructStatus(i);
834    if(status == CoinWarmStartBasis::atLowerBound) {
835      problem_->setxInit(i,colLow[i]);
836      problem_->setDualInit(i + numcols + numrows,0.);
837    }
838    else if(status == CoinWarmStartBasis::atUpperBound) {
839      problem_->setxInit(i,colUp[i]);
840      problem_->setDualInit(i + numrows,0.);
841    }
842    else {
843      problem_->setDualInit(i + numrows,0.);
844      problem_->setDualInit(i + numcols + numrows, 0.);
845    }
846  }
847  for(int i = 0; i < numrows ; i++) {
848    CoinWarmStartBasis::Status status = ws->getArtifStatus(i);
849    if(status == CoinWarmStartBasis::atLowerBound) {
850      problem_->setDualInit(i,0.);
851    }
852  }
853  int nElem = ws->values()->getNumElements();
854  const int * inds = ws->values()->getIndices();
855  const double * elems = ws->values()->getElements();
856
857  for(int i = 0 ; i < nElem ; i++) {
858    problem_->setxInit(inds[i],elems[i]);
859  }
860
861  if(IsValid(ws->warm_starter()))
862    problem_->SetWarmStarter(ws->warm_starter());
863  return 1;
864}
865
866static const char * OPT_SYMB="OPT";
867static const char * FAILED_SYMB="FAILED";
868static const char * INFEAS_SYMB="INFEAS";
869
870void
871IpoptInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
872    const char * whereFrom)
873{
874  totalNlpSolveTime_-=CoinCpuTime();
875  if(warmStarted)
876    optimization_status_ = app_->ReOptimizeTNLP(GetRawPtr(problem_));
877  else
878    optimization_status_ = app_->OptimizeTNLP(GetRawPtr(problem_));
879  totalNlpSolveTime_+=CoinCpuTime();
880  nCallOptimizeTNLP_++;
881  hasBeenOptimized_ = true;
882
883  //Options should have been printed if not done already turn off Ipopt output
884  if(!hasPrintedOptions) {
885    hasPrintedOptions = 1;
886    //app_->Options()->SetIntegerValue("print_level",0, true, true);
887    app_->Options()->SetStringValue("print_user_options","no", false, true);
888  }
889
890  const SmartPtr<SolveStatistics>  stats = app_->Statistics();
891
892#if 1
893  if(optimization_status_ == -10)//Too few degrees of freedom
894    {
895      std::cout<<"Too few degrees of freedom...."<<std::endl;
896      int numrows = getNumRows();
897      int numcols = getNumCols();
898
899      const double * colLower = getColLower();
900      const double * colUpper = getColUpper();
901
902
903      const double * rowLower = getRowLower();
904      const double * rowUpper = getRowUpper();
905
906      int numberFixed = 0;
907      for(int i = 0 ; i < numcols ; i++)
908        {
909          if(colUpper[i] - colLower[i] <= INT_BIAS)
910            {
911              numberFixed++;
912            }
913        }
914      int numberEqualities = 0;
915      for(int i = 0 ; i < numrows ; i++)
916        {
917          if(rowUpper[i] - rowLower[i] <= INT_BIAS)
918            {
919              numberEqualities++;
920            }     
921        }
922      if(numcols - numberFixed > numberEqualities)
923        {
924          throw UnsolvedError(optimization_status_);
925        }
926      double * saveColLow = CoinCopyOfArray(getColLower(), getNumCols());
927      double * saveColUp = CoinCopyOfArray(getColUpper(), getNumCols());
928
929      for(int i = 0; i < numcols && numcols - numberFixed <= numberEqualities ; i++)
930        {
931          if(colUpper[i] - colLower[i] <= INT_BIAS)
932            {
933              setColLower(i, saveColLow[i]-1e-06);
934              setColUpper(i, saveColLow[i]+1e-06);
935              numberFixed--;
936            }
937        }
938      solveAndCheckErrors(warmStarted, throwOnFailure, whereFrom);
939      //restore
940      for(int i = 0; i < numcols && numcols - numberFixed < numrows ; i++)
941        {
942          problem_->SetVariableLowerBound(i,saveColLow[i]);
943          problem_->SetVariableUpperBound(i,saveColUp[i]);
944        }
945      delete [] saveColLow;
946      delete [] saveColUp;
947      return;
948    }
949  else 
950#endif
951    if(optimization_status_ < -9)//Ipopt failed and the error can not be recovered, throw it
952  {
953    throw UnsolvedError(optimization_status_);
954  }
955
956  if(IsValid(stats)) {
957    totalIterations_ += stats->IterationCount();
958  }
959  else if (throwOnFailure)//something failed throw
960  {
961    throw SimpleError("No statistics available from Ipopt",whereFrom);
962  }
963  else return;
964
965  messageHandler()->message(IPOPT_SUMMARY, ipoptIMessages_)
966  <<whereFrom<<optimization_status_<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
967
968  if((nCallOptimizeTNLP_ % 20) == 1)
969    messageHandler()->message(LOG_HEAD, ipoptIMessages_)<<CoinMessageEol;
970
971
972  if (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) {
973    messageHandler()->message(SUSPECT_PROBLEM,
974        ipoptIMessages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
975    std::string subProbName;
976    getStrParam(OsiProbName, subProbName);
977    std::ostringstream os;
978    os<<"_"<<nCallOptimizeTNLP_;
979    subProbName+=os.str();
980    problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
981  }
982
983}
984
985////////////////////////////////////////////////////////////////////
986// Solve Methods                                                  //
987////////////////////////////////////////////////////////////////////
988/// Solve initial continuous relaxation
989void IpoptInterface::initialSolve()
990{
991  assert(IsValid(app_));
992  assert(IsValid(problem_));
993  if(!problem_->checkZeroDimension(optimization_status_)) {
994
995    if(!hasPrintedOptions) {
996      int printOptions;
997      app_->Options()->GetEnumValue("print_user_options",printOptions,"bonmin.");
998      if(printOptions)
999        app_->Options()->SetStringValue("print_user_options","yes");
1000    }
1001    solveAndCheckErrors(0,1,"initialSolve");
1002
1003    //Options should have been printed if not done already turn off Ipopt output
1004    if(!hasPrintedOptions) {
1005      hasPrintedOptions = 1;
1006      app_->Options()->SetStringValue("print_user_options","no");
1007      app_->Options()->SetIntegerValue("print_level",0);
1008    }
1009
1010    const char * status=OPT_SYMB;
1011    ;
1012    if(isAbandoned()) status=FAILED_SYMB;
1013    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1014    messageHandler()->message(LOG_FIRST_LINE, ipoptIMessages_)<<nCallOptimizeTNLP_
1015    <<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1016
1017    if(isAbandoned()) {
1018      resolveForRobustness(numRetryUnsolved_);
1019    }
1020    else if(numRetryInitial_)
1021    {
1022      resolveForCost(numRetryInitial_);
1023      /** Only do it once at the root.*/
1024      numRetryInitial_ = 0;
1025    }
1026  }
1027  else
1028    hasBeenOptimized_ = true;
1029}
1030
1031/** Resolve the continuous relaxation after problem modification.
1032 * \note for Ipopt, same as resolve */
1033void
1034IpoptInterface::resolve()
1035{
1036  if (!IsValid(app_->Statistics())) {
1037    messageHandler()->message(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, ipoptIMessages_)
1038    <<CoinMessageEol;
1039    initialSolve();
1040    return;
1041  }
1042  assert(IsValid(app_));
1043  assert(IsValid(problem_));
1044  if (INT_BIAS > 0.) {
1045    app_->Options()->SetStringValue("warm_start_same_structure", "yes");
1046  }
1047  else {
1048    app_->Options()->SetStringValue("warm_start_same_structure", "no");
1049  }
1050
1051  setWarmStartOptions();
1052
1053  if(!problem_->checkZeroDimension(optimization_status_)) {
1054    solveAndCheckErrors(1,1,"resolve");
1055
1056    const char * status=OPT_SYMB;
1057    ;
1058    if(isAbandoned()) status=FAILED_SYMB;
1059    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1060    messageHandler()->message(LOG_FIRST_LINE, ipoptIMessages_)<<nCallOptimizeTNLP_
1061    <<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1062
1063    if(isAbandoned()) {
1064      resolveForRobustness(numRetryUnsolved_);
1065    }
1066    else if(numRetryResolve_)
1067      resolveForCost(numRetryResolve_);
1068  }
1069  else
1070    hasBeenOptimized_ = true;
1071}
1072
1073void
1074IpoptInterface::resolveForCost(int numsolve)
1075{
1076
1077  /** Save current optimum. */
1078  double * point = new double[getNumCols()*3+ getNumRows()];
1079  double bestBound = getObjValue();
1080  CoinCopyN(getColSolution(),
1081      getNumCols(), point);
1082  CoinCopyN(getRowPrice(),
1083      2*getNumCols()+ getNumRows(),
1084      &point[getNumCols()]);
1085
1086  if(isProvenOptimal())
1087    messageHandler()->message(SOLUTION_FOUND,
1088        ipoptIMessages_)
1089    <<1<<getObjValue()<<bestBound
1090    <<CoinMessageEol;
1091  else
1092    messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
1093        ipoptIMessages_)
1094    <<1
1095    <<CoinMessageEol;
1096  for(int f = 0; f < numsolve ; f++) {
1097    messageHandler()->message(WARNING_RESOLVING,
1098        ipoptIMessages_)
1099    <<f+1<< CoinMessageEol ;
1100    randomStartingPoint();
1101    solveAndCheckErrors(0,0,"resolveForCost");
1102
1103
1104    const SmartPtr<SolveStatistics>  stats = app_->Statistics();
1105    if(IsValid(stats)) {
1106      totalIterations_ += stats->IterationCount();
1107    }
1108    else//something failed (random point generation is not that clever and this can happen (for ex asked to compute sqrt(-1))
1109      //just exit this procedure
1110    {
1111      return;
1112      //throw SimpleError("No statistics available from Ipopt","resolveForCost");
1113    }
1114
1115
1116
1117
1118
1119    const char * status=OPT_SYMB;
1120    ;
1121    char c=' ';
1122    if(isAbandoned()) {
1123      status=FAILED_SYMB;
1124    }
1125    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1126
1127
1128    //Is solution better than previous
1129    if(isProvenOptimal() &&
1130        getObjValue()<bestBound) {
1131      c='*';
1132      messageHandler()->message(BETTER_SOL, ipoptIMessages_)<<getObjValue()<<f+1<< CoinMessageEol;
1133      CoinCopyN(getColSolution(),
1134          getNumCols(), point);
1135      CoinCopyN(getRowPrice(),
1136          2*getNumCols()+ getNumRows(),
1137          &point[getNumCols()]);
1138      bestBound = getObjValue();
1139    }
1140
1141    messageHandler()->message(LOG_LINE, ipoptIMessages_)
1142    <<c<<f+1<<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1143
1144
1145    if(isProvenOptimal())
1146      messageHandler()->message(SOLUTION_FOUND,
1147          ipoptIMessages_)
1148      <<f+2<<getObjValue()<<bestBound
1149      <<CoinMessageEol;
1150    else if(!isAbandoned())
1151      messageHandler()->message(UNSOLVED_PROBLEM_FOUND,
1152          ipoptIMessages_)
1153      <<f+2
1154      <<CoinMessageEol;
1155    else
1156      messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
1157          ipoptIMessages_)
1158      <<f+2
1159      <<CoinMessageEol;
1160  }
1161  setColSolution(point);
1162  setRowPrice(&point[getNumCols()]);
1163  setWarmStartOptions();
1164  delete [] point;
1165
1166  optimization_status_ = app_->ReOptimizeTNLP(GetRawPtr(problem_));
1167  hasBeenOptimized_ = true;
1168}
1169
1170void
1171IpoptInterface::resolveForRobustness(int numsolve)
1172{
1173  //std::cerr<<"Resolving the problem for robustness"<<std::endl;
1174  //First remove warm start point and resolve
1175  unsetWarmStartOptions();
1176  messageHandler()->message(WARNING_RESOLVING,
1177      ipoptIMessages_)
1178  <<1<< CoinMessageEol ;
1179  solveAndCheckErrors(0,0,"resolveForRobustness");
1180
1181
1182  const char * status=OPT_SYMB;
1183  ;
1184  char c='*';
1185  if(isAbandoned()) {
1186    status=FAILED_SYMB;
1187    c=' ';
1188  }
1189  else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1190  messageHandler()->message(LOG_LINE, ipoptIMessages_)
1191  <<c<<1<<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1192
1193
1194  if(!isAbandoned()) {
1195    messageHandler()->message(WARN_SUCCESS_WS,
1196        ipoptIMessages_)
1197    << CoinMessageEol ;
1198    return; //we won go on
1199  }
1200
1201  //still unsolved try again with different random starting points
1202  for(int f = 0; f < numsolve ; f++) {
1203    messageHandler()->message(WARNING_RESOLVING,
1204        ipoptIMessages_)
1205    <<f+2<< CoinMessageEol ;
1206
1207    randomStartingPoint();
1208    solveAndCheckErrors(0,0,"resolveForRobustness");
1209
1210
1211    messageHandler()->message(IPOPT_SUMMARY, ipoptIMessages_)
1212    <<"resolveForRobustness"<<optimization_status_<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1213
1214
1215    const char * status=OPT_SYMB;
1216    ;
1217    char c='*';
1218    if(isAbandoned()) {
1219      status=FAILED_SYMB;
1220      c=' ';
1221    }
1222    else if(isProvenPrimalInfeasible()) status=INFEAS_SYMB;
1223    messageHandler()->message(LOG_LINE, ipoptIMessages_)
1224    <<c<<f+2<<status<<getObjValue()<<app_->Statistics()->IterationCount()<<app_->Statistics()->TotalCPUTime()<<CoinMessageEol;
1225
1226
1227    if(!isAbandoned()) {
1228      messageHandler()->message(WARN_SUCCESS_RANDOM,
1229          ipoptIMessages_)
1230      <<f+2
1231      << CoinMessageEol ;
1232      return; //we have found a solution and continue
1233    }
1234  }
1235  if(pretendFailIsInfeasible_) {
1236    if(pretendFailIsInfeasible_ == 1) {
1237      messageHandler()->message(WARN_CONTINUING_ON_FAILURE,
1238          ipoptIMessages_)
1239      <<CoinMessageEol;
1240      hasContinuedAfterNlpFailure_ = 1;
1241    }
1242    return;
1243  }
1244  else {
1245    throw UnsolvedError(optimization_status_);
1246  }
1247}
1248
1249////////////////////////////////////////////////////////////////////
1250// Methods returning info on how the solution process terminated  //
1251////////////////////////////////////////////////////////////////////
1252/// Are there a numerical difficulties?
1253bool IpoptInterface::isAbandoned() const
1254{
1255  return (
1256        (optimization_status_==Ipopt::Maximum_Iterations_Exceeded)||
1257        (optimization_status_==Ipopt::Restoration_Failed)||
1258        (optimization_status_==Ipopt::Error_In_Step_Computation)||
1259        (optimization_status_==Ipopt::Not_Enough_Degrees_Of_Freedom)||
1260        (optimization_status_==Ipopt::Invalid_Problem_Definition)||
1261        (optimization_status_==Ipopt::Invalid_Option)||
1262        (optimization_status_==Ipopt::Invalid_Number_Detected)||
1263        (optimization_status_==Ipopt::Unrecoverable_Exception)||
1264        (optimization_status_==Ipopt::NonIpopt_Exception_Thrown)||
1265        (optimization_status_==Ipopt::Insufficient_Memory)||
1266        (optimization_status_==Ipopt::Internal_Error)
1267      );
1268}
1269
1270/// Is optimality proven?
1271bool IpoptInterface::isProvenOptimal() const
1272{
1273  if (optimization_status_==Ipopt::Search_Direction_Becomes_Too_Small) {
1274    std::cerr<<"Warning : need to verify that Search_Direction_Becomes_Too_Small is indeed Ok"<<std::endl;
1275  }
1276  return (optimization_status_==Ipopt::Solve_Succeeded ||
1277      optimization_status_==Ipopt::Search_Direction_Becomes_Too_Small ||
1278      optimization_status_==Ipopt::Solved_To_Acceptable_Level);
1279}
1280/// Is primal infeasiblity proven?
1281bool IpoptInterface::isProvenPrimalInfeasible() const
1282{
1283  return (optimization_status_ == Ipopt::Infeasible_Problem_Detected);
1284}
1285/// Is dual infeasiblity proven?
1286bool IpoptInterface::isProvenDualInfeasible() const
1287{
1288  throw SimpleError("Don't have this optimization status yet.",
1289      "isProvenDualInfeasible");
1290}
1291/// Is the given primal objective limit reached?
1292bool IpoptInterface::isPrimalObjectiveLimitReached() const
1293{
1294  std::cerr<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<std::endl;
1295  return 0;
1296}
1297/// Is the given dual objective limit reached?
1298bool IpoptInterface::isDualObjectiveLimitReached() const
1299{
1300  //  std::cerr<<"Warning : isDualObjectiveLimitReached not implemented yet"<<std::endl;
1301  return (optimization_status_==Ipopt::Diverging_Iterates);
1302
1303}
1304/// Iteration limit reached?
1305bool IpoptInterface::isIterationLimitReached() const
1306{
1307  return (optimization_status_==Ipopt::Maximum_Iterations_Exceeded);
1308  //  return (problem_->optimization_status()==Ipopt::Maximum_Iterations_Exceeded);
1309}
1310////////////////////////////////////////////////////////////////////
1311// Problem information methods                                    //
1312////////////////////////////////////////////////////////////////////
1313/// Get number of columns
1314int IpoptInterface::getNumCols() const
1315{
1316
1317  return problem_->num_variables();
1318}
1319
1320
1321/// Get number of rows
1322int
1323IpoptInterface::getNumRows() const
1324{
1325  return problem_->num_constraints();
1326}
1327
1328const double *
1329IpoptInterface::getColLower() const
1330{
1331  return problem_->x_l();
1332}
1333
1334const double *
1335IpoptInterface::getColUpper() const
1336{
1337  return problem_->x_u();
1338}
1339
1340void
1341IpoptInterface::readVarNames() const
1342{
1343  delete []varNames_;
1344  varNames_ = NULL;
1345  std::string probName;
1346  getStrParam(OsiProbName, probName);
1347  IpCbcColReader colRead(probName);
1348  if(colRead.readFile()) {
1349    varNames_ = new std::string[problem_->num_variables()];
1350    colRead.copyNames(varNames_,problem_->num_variables());
1351    hasVarNamesFile_ = true;
1352  }
1353  else
1354    hasVarNamesFile_ = false;
1355}
1356
1357///get name of a variable
1358const std::string *
1359IpoptInterface::getVarNames() const
1360{
1361  if(varNames_ == NULL && hasVarNamesFile_ ) {
1362    readVarNames();
1363  }
1364  return varNames_;
1365}
1366
1367void IpoptInterface::extractSenseRhsAndRange() const
1368{
1369  assert(rowsense_==NULL&&rhs_==NULL&&rowrange_==NULL);
1370  int numrows = problem_->num_constraints();
1371  const double * rowLower = getRowLower();
1372  const double * rowUpper = getRowUpper();
1373  rowsense_ = new char [numrows];
1374  rhs_ = new double [numrows];
1375  rowrange_ = new double [numrows];
1376  for(int i = 0 ; i < numrows ; i++) {
1377    rowrange_[i]=0.;
1378    convertBoundToSense(rowLower[i], rowUpper[i], rowsense_[i], rhs_[i], rowrange_[i]);
1379  }
1380}
1381
1382/** Get pointer to array[getNumRows()] of row constraint senses.
1383    <ul>
1384    <li>'L': <= constraint
1385    <li>'E': =  constraint
1386    <li>'G': >= constraint
1387    <li>'R': ranged constraint
1388    <li>'N': free constraint
1389    </ul>
1390*/
1391const char *
1392IpoptInterface::getRowSense() const
1393{
1394  if(rowsense_==NULL) {
1395    extractSenseRhsAndRange();
1396  }
1397  return rowsense_;
1398}
1399
1400/** Get pointer to array[getNumRows()] of rows right-hand sides
1401    <ul>
1402    <li> if rowsense()[i] == 'L' then rhs()[i] == rowupper()[i]
1403    <li> if rowsense()[i] == 'G' then rhs()[i] == rowlower()[i]
1404    <li> if rowsense()[i] == 'R' then rhs()[i] == rowupper()[i]
1405    <li> if rowsense()[i] == 'N' then rhs()[i] == 0.0
1406    </ul>
1407*/
1408const double *
1409IpoptInterface::getRightHandSide() const
1410{
1411  if(rhs_==NULL) {
1412    extractSenseRhsAndRange();
1413  }
1414  return rhs_;
1415}
1416
1417/** Get pointer to array[getNumRows()] of row ranges.
1418    <ul>
1419    <li> if rowsense()[i] == 'R' then
1420    rowrange()[i] == rowupper()[i] - rowlower()[i]
1421    <li> if rowsense()[i] != 'R' then
1422    rowrange()[i] is 0.0
1423    </ul>
1424*/
1425const double *
1426IpoptInterface::getRowRange() const
1427{
1428  if(rowrange_==NULL) {
1429    extractSenseRhsAndRange();
1430  }
1431  return rowrange_;
1432}
1433
1434const double *
1435IpoptInterface::getRowLower() const
1436{
1437  return problem_->g_l();
1438}
1439
1440const double *
1441IpoptInterface::getRowUpper() const
1442{
1443  return problem_->g_u();
1444}
1445
1446/// Return true if column is continuous
1447bool
1448IpoptInterface::isContinuous(int colNumber) const
1449{
1450  return (problem_->var_types()[colNumber]==Ipopt::TMINLP::CONTINUOUS);
1451}
1452
1453/// Return true if column is binary
1454bool
1455IpoptInterface::isBinary(int colNumber) const
1456{
1457  return (problem_->var_types()[colNumber]==Ipopt::TMINLP::BINARY);
1458}
1459
1460/** Return true if column is integer.
1461    Note: This function returns true if the the column
1462    is binary or a general integer.
1463*/
1464bool
1465IpoptInterface::isInteger(int colNumber) const
1466{
1467  return ((problem_->var_types()[colNumber]==Ipopt::TMINLP::BINARY)||
1468      (problem_->var_types()[colNumber]==Ipopt::TMINLP::INTEGER));
1469}
1470
1471/// Return true if column is general integer
1472bool
1473IpoptInterface::isIntegerNonBinary(int colNumber) const
1474{
1475  return (problem_->var_types()[colNumber]==Ipopt::TMINLP::INTEGER);
1476}
1477/// Return true if column is binary and not fixed at either bound
1478bool
1479IpoptInterface::isFreeBinary(int colNumber) const
1480{
1481  return ((problem_->var_types()[colNumber]==Ipopt::TMINLP::BINARY)
1482      &&((getColUpper()[colNumber]-getColLower()[colNumber]) > 1 - 1e-09));
1483}
1484
1485/// Get solver's value for infinity
1486double
1487IpoptInterface::getInfinity() const
1488{
1489  return 1e100;
1490}
1491
1492/// Get pointer to array[getNumCols()] of primal solution vector
1493const double *
1494IpoptInterface::getColSolution() const
1495{
1496  if(hasBeenOptimized_)
1497    return problem_->x_sol();
1498  else
1499    return problem_->x_init();
1500}
1501
1502/// Get pointer to array[getNumRows()] of dual prices
1503const double *
1504IpoptInterface::getRowPrice() const
1505{
1506  if(hasBeenOptimized_)
1507    return problem_->duals_sol();
1508  else
1509    return problem_->duals_init();
1510}
1511
1512/// Get a pointer to array[getNumCols()] of reduced costs
1513const double *
1514IpoptInterface::getReducedCost() const
1515{
1516  std::cerr<<"WARNING : trying to access reduced cost in Ipopt always retrun 0"<<std::endl;
1517  if(reducedCosts_==NULL) {
1518    reducedCosts_ = new double [getNumCols()];
1519    CoinFillN(reducedCosts_,getNumCols(),0.);
1520  }
1521  return reducedCosts_;
1522}
1523
1524/** Get pointer to array[getNumRows()] of row activity levels (constraint
1525    matrix times the solution vector */
1526const double *
1527IpoptInterface::getRowActivity() const
1528{
1529  return problem_->g_sol();
1530}
1531
1532/** Get how many iterations it took to solve the problem (whatever
1533    "iteration" mean to the solver.
1534*/
1535int
1536IpoptInterface::getIterationCount() const
1537{
1538  if(IsValid(app_->Statistics()))
1539    return app_->Statistics()->IterationCount();
1540  else
1541    return 0;
1542}
1543
1544
1545/** Set a single column lower bound.
1546    Use -getInfinity() for -infinity. */
1547void
1548IpoptInterface::setColLower( int elementIndex, double elementValue )
1549{
1550  //  if(fabs(problem_->x_l()[elementIndex]-elementValue)>1e-06)
1551  problem_->SetVariableLowerBound(elementIndex,elementValue);
1552  hasBeenOptimized_ = false;
1553}
1554
1555/** Set a single column upper bound.
1556    Use getInfinity() for infinity. */
1557void
1558IpoptInterface::setColUpper( int elementIndex, double elementValue )
1559{
1560  //  if(fabs(problem_->x_u()[elementIndex]-elementValue)>1e-06)
1561  problem_->SetVariableUpperBound(elementIndex,elementValue);
1562  hasBeenOptimized_ = false;
1563}
1564
1565/** Set a single row lower bound.
1566    Use -getInfinity() for -infinity. */
1567void
1568IpoptInterface::setRowLower( int elementIndex, double elementValue )
1569{
1570  throw SimpleError("Not implemented yet but should be if necessary.",
1571      "setRowLower");
1572  hasBeenOptimized_ = false;
1573}
1574
1575/** Set a single row upper bound.
1576    Use getInfinity() for infinity. */
1577void
1578IpoptInterface::setRowUpper( int elementIndex, double elementValue )
1579{
1580  throw SimpleError("Not implemented yet but should be if necessary.",
1581      "setRowUpper");
1582  hasBeenOptimized_ = false;
1583}
1584
1585/** Set the type of a single row */
1586void
1587IpoptInterface::setRowType(int index, char sense, double rightHandSide,
1588    double range)
1589{
1590  throw SimpleError("Not implemented yet but should be if necessary.",
1591      "setRowType");
1592  hasBeenOptimized_ = false;
1593}
1594
1595
1596/// Set the objective function sense.
1597/// (1 for min (default), -1 for max)
1598void
1599IpoptInterface::setObjSense(double s)
1600{
1601  throw SimpleError("Can not change objective sense of an Ipopt problem.",
1602      "setObjSense");
1603  hasBeenOptimized_ = false;
1604}
1605
1606/** Set the primal solution variable values
1607
1608colsol[getNumCols()] is an array of values for the primal variables.
1609These values are copied to memory owned by the solver interface object
1610or the solver.  They will be returned as the result of getColSolution()
1611until changed by another call to setColSolution() or by a call to any
1612solver routine.  Whether the solver makes use of the solution in any
1613way is solver-dependent.
1614*/
1615void
1616IpoptInterface::setColSolution(const double *colsol)
1617{
1618  problem_->setxInit(getNumCols(), colsol);
1619  hasBeenOptimized_ = false;
1620  //  problem_->SetStartingPoint(getNumCols(), colsol);
1621}
1622
1623/** Set dual solution variable values
1624
1625rowprice[getNumRows()] is an array of values for the dual
1626variables. These values are copied to memory owned by the solver
1627interface object or the solver.  They will be returned as the result of
1628getRowPrice() until changed by another call to setRowPrice() or by a
1629call to any solver routine.  Whether the solver makes use of the
1630solution in any way is solver-dependent.
1631*/
1632
1633void
1634IpoptInterface::setRowPrice(const double * rowprice)
1635{
1636  problem_->setDualsInit(getNumCols()*2 + getNumRows(), rowprice);
1637  hasBeenOptimized_ = false;
1638}
1639
1640
1641/** Set the index-th variable to be a continuous variable */
1642void
1643IpoptInterface::setContinuous(int index)
1644{
1645  problem_->SetVariableType(index, Ipopt::TMINLP::CONTINUOUS);
1646  hasBeenOptimized_ = false;
1647}
1648/** Set the index-th variable to be an integer variable */
1649void
1650IpoptInterface::setInteger(int index)
1651{
1652  problem_->SetVariableType(index, Ipopt::TMINLP::INTEGER);
1653  hasBeenOptimized_ = false;
1654}
1655
1656/// Get objective function value (can't use default)
1657double
1658IpoptInterface::getObjValue() const
1659{
1660  return problem_->obj_value();
1661}
1662
1663
1664CoinWarmStart* IpoptInterface::getEmptyWarmStart () const
1665{
1666  return (dynamic_cast<CoinWarmStart *>(new IpoptWarmStart(1)));
1667}
1668
1669//#############################################################################
1670// Parameter related methods
1671//#############################################################################
1672
1673bool
1674IpoptInterface::setIntParam(OsiIntParam key, int value)
1675{
1676  //  debugMessage("OsiCpxSolverInterface::setIntParam(%d, %d)\n", key, value);
1677
1678  bool retval = false;
1679  switch (key) {
1680  case OsiMaxNumIteration:
1681    retval = false;
1682    break;
1683  case OsiMaxNumIterationHotStart:
1684    if( value >= 0 ) {
1685      retval = false;
1686    }
1687    else
1688      retval = false;
1689    break;
1690  case OsiLastIntParam:
1691    retval = false;
1692    break;
1693  }
1694  return retval;
1695}
1696
1697//-----------------------------------------------------------------------------
1698
1699bool
1700IpoptInterface::setDblParam(OsiDblParam key, double value)
1701{
1702  //  debugMessage("IpoptInterface::setDblParam(%d, %g)\n", key, value);
1703
1704  bool retval = false;
1705  switch (key) {
1706  case OsiDualObjectiveLimit:
1707    OsiDualObjectiveLimit_ = value;
1708    retval = true;
1709    break;
1710  case OsiPrimalObjectiveLimit:
1711    std::cerr<<"Can not set primal objective limit parameter"<<std::endl;
1712    retval = false;
1713    break;
1714  case OsiDualTolerance:
1715    std::cerr<<"Can not set dual tolerance parameter"<<std::endl;
1716    retval = false;
1717    break;
1718  case OsiPrimalTolerance:
1719    std::cerr<<"Can not set primal tolerance parameter"<<std::endl;
1720    retval = false;
1721  case OsiObjOffset:
1722    retval = OsiSolverInterface::setDblParam(key,value);
1723    break;
1724  case OsiLastDblParam:
1725    retval = false;
1726    break;
1727  }
1728  return retval;
1729}
1730
1731
1732//-----------------------------------------------------------------------------
1733
1734bool
1735IpoptInterface::setStrParam(OsiStrParam key, const std::string & value)
1736{
1737  //  debugMessage("IpoptInterface::setStrParam(%d, %s)\n", key, value.c_str());
1738
1739  bool retval=false;
1740  switch (key) {
1741  case OsiProbName:
1742    OsiSolverInterface::setStrParam(key,value);
1743    return retval = true;
1744  case OsiSolverName:
1745    return false;
1746  case OsiLastStrParam:
1747    return false;
1748  }
1749  return false;
1750}
1751
1752//-----------------------------------------------------------------------------
1753
1754bool
1755IpoptInterface::getIntParam(OsiIntParam key, int& value) const
1756{
1757  //  debugMessage("IpoptInterface::getIntParam(%d)\n", key);
1758
1759  bool retval = false;
1760  switch (key) {
1761  case OsiMaxNumIteration:
1762    retval = false;
1763    break;
1764  case OsiMaxNumIterationHotStart:
1765    retval = false;
1766    break;
1767  case OsiLastIntParam:
1768    retval = false;
1769    break;
1770  }
1771  return retval;
1772}
1773
1774//-----------------------------------------------------------------------------
1775
1776bool
1777IpoptInterface::getDblParam(OsiDblParam key, double& value) const
1778{
1779  //  debugMessage("IpoptInterface::getDblParam(%d)\n", key);
1780
1781  bool retval = false;
1782  switch (key) {
1783  case OsiDualObjectiveLimit:
1784    value = OsiDualObjectiveLimit_;
1785    retval = true;
1786    break;
1787  case OsiPrimalObjectiveLimit:
1788    value = getInfinity();
1789    retval = true;
1790    break;
1791  case OsiDualTolerance:
1792    retval = false;
1793    break;
1794  case OsiPrimalTolerance:
1795    retval = false;
1796    break;
1797  case OsiObjOffset:
1798    retval = OsiSolverInterface::getDblParam(key, value);
1799    break;
1800  case OsiLastDblParam:
1801    retval = false;
1802    break;
1803  }
1804  return retval;
1805}
1806
1807
1808//-----------------------------------------------------------------------------
1809
1810bool
1811IpoptInterface::getStrParam(OsiStrParam key, std::string & value) const
1812{
1813  //  debugMessage("IpoptInterface::getStrParam(%d)\n", key);
1814
1815  switch (key) {
1816  case OsiProbName:
1817    OsiSolverInterface::getStrParam(key, value);
1818    break;
1819  case OsiSolverName:
1820    value = "Ipopt";
1821    break;
1822  case OsiLastStrParam:
1823    return false;
1824  }
1825
1826  return true;
1827}
1828
1829void
1830IpoptInterface::turnOffIpoptOutput()
1831{
1832  std::string opt="print_level";
1833  app_->Options()->SetIntegerValue(opt, (Index)J_NONE,true);
1834}
1835void
1836IpoptInterface::turnOnIpoptOutput()
1837{
1838  std::string opt="print_level";
1839  app_->Options()->SetIntegerValue(opt, (Index)J_SUMMARY,true);
1840}
1841
1842void
1843IpoptInterface::randomStartingPoint()
1844{
1845  int numcols = getNumCols();
1846  const double * colLower = getColLower();
1847  const double * colUpper = getColUpper();
1848  double * sol = new double[numcols];
1849  for(int i = 0 ; i < numcols ; i++) {
1850    double lower = min(-maxRandomRadius_,colUpper[i] - maxRandomRadius_);
1851    lower = max(colLower[i], lower);
1852    double upper = max(maxRandomRadius_,colLower[i] + maxRandomRadius_);
1853    upper = min(colUpper[i],upper);
1854    lower = min(upper,lower);
1855    upper = max(upper, lower);
1856    double interval = upper - lower;
1857    sol[i] = CoinDrand48()*(interval) + lower;
1858//    printf("%f in [%f,%f]\n",sol[i],lower,upper);
1859    //  std::cout<<interval<<"\t";
1860  }
1861  //std::cout<<std::endl;
1862  unsetWarmStartOptions();
1863  setColSolution(sol);
1864  delete [] sol;
1865}
1866
1867/*******************************************************************************/
1868// Class for throwing errors reported from Ipopt
1869/******************************************************************************/
1870
1871std::string
1872IpoptInterface::UnsolvedError::errorNames[17] ={"Solve succeeded",
1873    "Solved to acceptable level",
1874    "Infeasible Problem Detected",
1875    "Search_Direction_Becomes_Too_Small",
1876    "Diverging iterates",
1877    "User Requested Stop",
1878    "Maximum Iterations Exceeded",
1879    "Restoration Failed",
1880    "Error in step computation",
1881    "Not enough degrees of freedom",
1882    "Invalid problem definition",
1883    "Invalid option",
1884    "Invalid number detected",
1885    "Unrecoverable exception",
1886    "NonIpopt exception thrown",
1887    "Insufficient memory",
1888    "Internal error"};
1889
1890IpoptInterface::UnsolvedError::UnsolvedError(int errorNum)
1891{
1892  errorNum_ = errorNum;
1893}
1894
1895const std::string &
1896IpoptInterface::UnsolvedError::errorName() const
1897{
1898  if(errorNum_ >=0)
1899    return errorNames[errorNum_];
1900  if(errorNum_ == -1) return errorNames[6];
1901  else if(errorNum_ == -2) return errorNames[7];
1902  else if(errorNum_ == -3) return errorNames[8];
1903  else if(errorNum_ == -10) return errorNames[9];
1904  else if(errorNum_ == -11) return errorNames[10];
1905  else if(errorNum_ == -12) return errorNames[11];
1906  else if(errorNum_ == -13) return errorNames[12];
1907  else if(errorNum_ == -100) return errorNames[13];
1908  else if(errorNum_ == -101) return errorNames[14];
1909  else if(errorNum_ == -102) return errorNames[15];
1910  else if(errorNum_ == -199) return errorNames[16];
1911  throw IpoptInterface::SimpleError("UnsolvedError::errorName()","Unrecognized optimization status in ipopt.");
1912}
1913
1914void
1915IpoptInterface::UnsolvedError::printError(std::ostream &os)
1916{
1917  os<<"Ipopt exited with error code "<<errorNum_<<" "<<errorName()<<std::endl;
1918}
1919
1920
1921/** This methods initialiaze arrays for storing the jacobian */
1922int IpoptInterface::initializeJacobianArrays()
1923{
1924  Ipopt::Index n, m, nnz_h_lag;
1925  Ipopt::TNLP::IndexStyleEnum index_style;
1926  tminlp_->get_nlp_info( n, m, nnz_jac, nnz_h_lag, index_style);
1927
1928  if(jRow_ != NULL) delete jRow_;
1929  if(jCol_ != NULL) delete jCol_;
1930  if(jValues_ != NULL) delete jValues_;
1931
1932  jRow_ = new Ipopt::Index[nnz_jac];
1933  jCol_ = new Ipopt::Index[nnz_jac];
1934  jValues_ = new Ipopt::Number[nnz_jac];
1935  tminlp_->eval_jac_g(n, NULL, 0, m, nnz_jac, jRow_, jCol_, NULL);
1936
1937
1938  if(constTypes_ != NULL) delete [] constTypes_;
1939  if(constTypesNum_ != NULL) delete [] constTypesNum_;
1940
1941  constTypes_ = new Ipopt::TMINLP::ConstraintType[getNumRows()];
1942  tminlp_->get_constraints_types(getNumRows(), constTypes_);
1943  constTypesNum_ = new int[getNumRows()];
1944  for(int i = 0; i < getNumRows() ; i++) {
1945    if(constTypes_[i]==Ipopt::TMINLP::LINEAR) {
1946      constTypesNum_[i] = nLinear_++;
1947    }
1948    else if(constTypes_[i]==Ipopt::TMINLP::NON_LINEAR) {
1949      constTypesNum_[i] = nNonLinear_++;
1950    }
1951  }
1952  return nnz_jac;
1953}
1954
1955
1956/** get NLP constraint violation of current point */
1957double
1958IpoptInterface::getConstraintViolation()
1959{
1960  int numcols = getNumCols();
1961  int numrows = getNumRows();
1962  double * g = new double[numrows];
1963  tminlp_->eval_g(numcols,getColSolution(),1,numrows,g);
1964  const double * rowLower = getRowLower();
1965  const double * rowUpper = getRowUpper();
1966
1967
1968  double norm = 0;
1969  for(int i = 0; i< numrows ; i++) {
1970    if(constTypes_[i] == Ipopt::TMINLP::NON_LINEAR) {
1971      double rowViolation = max(0.,rowLower[i] - g[i]);
1972      if(rowViolation  > 0.) rowViolation  /= fabs(rowLower[i]);
1973      double rowViolation2 = max(0.,g[i] - rowUpper[i]);
1974      if(rowViolation2 > 0.) rowViolation2 /= fabs(rowUpper[i]);
1975      norm = max(rowViolation, norm);
1976    }
1977  }
1978  return norm;
1979}
1980
1981
1982//A procedure to try to remove small coefficients in OA cuts (or make it non small
1983static inline
1984bool cleanNnz(double &value, double colLower, double colUpper,
1985    double rowLower, double rowUpper, double colsol,
1986    double & lb, double &ub, double tiny, double veryTiny)
1987{
1988  if(fabs(value)>= tiny) return 1;
1989
1990  if(fabs(value)<veryTiny) return 0;//Take the risk?
1991
1992  //try and remove
1993  double infty = 1e20;
1994  bool colUpBounded = colUpper < 10000;
1995  bool colLoBounded = colLower > -10000;
1996  bool rowNotLoBounded =  rowLower <= - infty;
1997  bool rowNotUpBounded = rowUpper >= infty;
1998  bool pos =  value > 0;
1999
2000  if(!rowNotLoBounded && ! rowNotUpBounded)//would have to either choose side or duplicate cut
2001  {
2002    std::cerr<<"OA on non-convex constraint is very experimental, don't know how to remove"
2003    <<std::endl;
2004  }
2005
2006  if(colLoBounded && pos && rowNotUpBounded) {
2007    lb += value * (colsol - colLower);
2008    return 0;
2009  }
2010  else
2011    if(colLoBounded && !pos && rowNotLoBounded) {
2012      ub += value * (colsol - colLower);
2013      return 0;
2014    }
2015    else
2016      if(colUpBounded && !pos && rowNotUpBounded) {
2017        lb += value * (colsol - colUpper);
2018        return 0;
2019      }
2020      else
2021        if(colUpBounded && pos && rowNotLoBounded) {
2022          ub += value * (colsol - colUpper);
2023          return 0;
2024        }
2025  //can not remove coefficient increase it to smallest non zero
2026  if(pos) value = tiny;
2027  else
2028    value = - tiny;
2029  return 1;
2030}
2031
2032/** Get the outer approximation constraints at the current optimum of the
2033  current ipopt problem */
2034void
2035IpoptInterface::getOuterApproximation(OsiCuts &cs, bool getObj)
2036{
2037  int n,m, nnz_jac_g, nnz_h_lag;
2038  Ipopt::TNLP::IndexStyleEnum index_style;
2039  tminlp_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2040  if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2041    initializeJacobianArrays();
2042  assert(jRow_ != NULL);
2043  assert(jCol_ != NULL);
2044  double * g = new double[m];
2045  tminlp_->eval_jac_g(n, getColSolution(), 1, m, nnz_jac_g, NULL, NULL, jValues_);
2046  tminlp_->eval_g(n,getColSolution(),1,m,g);
2047  //As jacobian is stored by cols fill OsiCuts with cuts
2048  CoinPackedVector * cuts = new CoinPackedVector[nNonLinear_ + 1];
2049  double * lb = new double[nNonLinear_ + 1];
2050  double * ub = new double[nNonLinear_ + 1];
2051  const double * rowLower = getRowLower();
2052  const double * rowUpper = getRowUpper();
2053  const double * colLower = getColLower();
2054  const double * colUpper = getColUpper();
2055
2056  for(int i = 0; i< m ; i++) {
2057    if(constTypes_[i] == Ipopt::TMINLP::NON_LINEAR) {
2058      lb[constTypesNum_[i]] = rowLower[i] - g[i];
2059      ub[constTypesNum_[i]] = rowUpper[i] - g[i];
2060    }
2061  }
2062//double infty = getInfinity();
2063  for(int i = 0 ; i < nnz_jac_g ; i++) {
2064    if(constTypes_[jRow_[i] - 1] == Ipopt::TMINLP::NON_LINEAR) {
2065      //"clean" coefficient
2066      if(cleanNnz(jValues_[i],colLower[jCol_[i] - 1], colUpper[jCol_[i]-1],
2067          rowLower[jRow_[i] - 1], rowUpper[jRow_[i] - 1],
2068          getColSolution()[jCol_[i] - 1],
2069          lb[constTypesNum_[jRow_[i] - 1]],
2070          ub[constTypesNum_[jRow_[i] - 1]], tiny_, veryTiny_)) {
2071        cuts[constTypesNum_[jRow_[i] - 1]].insert(jCol_[i]-1,jValues_[i]);
2072        lb[constTypesNum_[jRow_[i] - 1]] += jValues_[i] * getColSolution()[jCol_ [i] - 1];
2073        ub[constTypesNum_[jRow_[i] - 1]] += jValues_[i] * getColSolution()[jCol_ [i] - 1];
2074      }
2075    }
2076  }
2077
2078  for(int i = 0; i< nNonLinear_ ; i++) {
2079    OsiRowCut newCut;
2080    //    if(lb[i]>-1e20) assert (ub[i]>1e20);
2081
2082    newCut.setGloballyValid();
2083    newCut.setLb(lb[i]);
2084    newCut.setUb(ub[i]);
2085    newCut.setRow(cuts[i]);
2086    //    CftValidator validator;
2087    //    validator(newCut);
2088    cs.insert(newCut);
2089  }
2090
2091  delete[] g;
2092  delete [] cuts;
2093
2094  if(getObj)  // Get the objective cuts
2095  {
2096    double * obj = new double [n];
2097    tminlp_->eval_grad_f(n, getColSolution(), 1,obj);
2098    double f;
2099    tminlp_->eval_f(n, getColSolution(), 1, f);
2100
2101    CoinPackedVector v;
2102    v.reserve(n);
2103    lb[nNonLinear_] = -f;
2104    ub[nNonLinear_] = -f;
2105    //double minCoeff = 1e50;
2106    for(int i = 0; i<n ; i++)
2107    {
2108      if(cleanNnz(obj[i],colLower[i], colUpper[i],
2109          -getInfinity(), 0,
2110          getColSolution()[i],
2111          lb[nNonLinear_],
2112          ub[nNonLinear_],tiny_, 1e-15)) {
2113        //            minCoeff = min(fabs(obj[i]), minCoeff);
2114        v.insert(i,obj[i]);
2115        lb[nNonLinear_] += obj[i] * getColSolution()[i];
2116        ub[nNonLinear_] += obj[i] * getColSolution()[i];
2117      }
2118    }
2119    v.insert(n,-1);
2120    OsiRowCut newCut;
2121    newCut.setGloballyValid();
2122    newCut.setRow(v);
2123    newCut.setLb(-DBL_MAX/*Infinity*/);
2124    newCut.setUb(ub[nNonLinear_]);
2125//     CftValidator validator;
2126//     validator(newCut);
2127    cs.insert(newCut);
2128    delete [] obj;
2129  }
2130
2131  // setColSolution(problem()->x_sol());
2132  //  setRowPrice(problem()->duals_sol());
2133  delete []lb;
2134  delete[]ub;
2135}
2136
2137double
2138IpoptInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs)
2139{
2140  if(! IsValid(feasibilityProblem_)) {
2141    throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2142  }
2143  feasibilityProblem_->set_dist2point_obj(n,(const Ipopt::Number *) x_bar,(const Ipopt::Index *) inds);
2144  nCallOptimizeTNLP_++;
2145  totalNlpSolveTime_-=CoinCpuTime();
2146  Ipopt::IpoptApplication app2;
2147  app2.Options()->SetIntegerValue("print_level", (Ipopt::Index) 0);
2148  optimization_status_ = app2.OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2149  totalNlpSolveTime_+=CoinCpuTime();
2150  getOuterApproximation(cs, 0);
2151  hasBeenOptimized_=true;
2152  return getObjValue();
2153}
2154
2155
2156
2157
2158
2159void
2160IpoptInterface::extractLinearRelaxation(OsiSolverInterface &si, bool getObj)
2161{
2162  initialSolve();
2163
2164  CoinPackedMatrix mat;
2165  double * rowLow = NULL;
2166  double * rowUp = NULL;
2167
2168  int n;
2169  int m;
2170  int nnz_jac_g;
2171  int nnz_h_lag;
2172  Ipopt::TNLP::IndexStyleEnum index_style;
2173  //Get problem information
2174  tminlp_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2175
2176  //if not allocated allocate spaced for stroring jacobian
2177  if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2178    initializeJacobianArrays();
2179
2180  //get Jacobian
2181  tminlp_->eval_jac_g(n, getColSolution(), 1, m, nnz_jac_g, NULL, NULL, jValues_);
2182
2183
2184  double *g = new double[m];
2185  tminlp_->eval_g(n, getColSolution(),1,m,g);
2186
2187  rowLow = new double[m];
2188  rowUp = new double[m];
2189  const double * rowLower = getRowLower();
2190  const double * rowUpper = getRowUpper();
2191  const double * colLower = getColLower();
2192  const double * colUpper = getColUpper();
2193  assert(m==getNumRows());
2194  for(int i = 0 ; i < m ; i++) {
2195    {
2196      if(constTypes_[i] == Ipopt::TMINLP::NON_LINEAR) {
2197        rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
2198        rowUp[i] =  (rowUpper[i] - g[i]) + 1e-07;
2199      }
2200      else {
2201        rowLow[i] = (rowLower[i] - g[i]);
2202        rowUp[i] =  (rowUpper[i] - g[i]);
2203      }
2204    }
2205  }
2206
2207  //Then convert everything to a CoinPackedMatrix (col ordered)
2208  CoinBigIndex * inds = new CoinBigIndex[nnz_jac_g];
2209  double * vals = new double [nnz_jac_g];
2210  int * start = new int[n+1];
2211  int * length = new int[n];
2212  bool needOrder = false;
2213  int nnz = 0;
2214  for(int k = 0; k < jCol_[0]; k++) {
2215    start[k] = nnz;
2216    length[k] = 0;
2217  }
2218  for(int i = 0 ; i < nnz_jac_g - 1 ; i++) {
2219    {
2220      if(jCol_[i + 1] < jCol_ [i] || ( jCol_ [i+1] == jCol_[i] && jRow_[i + 1] <= jRow_[i]) ) {
2221        needOrder = 1;
2222        break;
2223      }
2224      if(Ipopt::TMINLP::LINEAR //Always accept coefficients from linear constraints
2225          || //For other clean tinys
2226          cleanNnz(jValues_[i],colLower[jCol_[i] - 1], colUpper[jCol_[i]-1],
2227              rowLower[jRow_[i] - 1], rowUpper[jRow_[i] - 1],
2228              getColSolution()[jCol_[i] - 1],
2229              rowLow[jRow_[i] - 1],
2230              rowUp[jRow_[i] -1], tiny_, veryTiny_)) {
2231        vals[nnz] = jValues_[i];
2232        rowLow[jRow_[i] - 1] += jValues_[i] * getColSolution()[jCol_ [i] - 1];
2233        rowUp[jRow_[i] -1] += jValues_[i] *getColSolution()[jCol_[i] -1];
2234        inds[nnz] = jRow_[i ] - 1;
2235        length[jCol_[i] - 1]++;
2236        nnz++;
2237      }
2238    }
2239    if(jCol_[i + 1] > jCol_[i]) {
2240      for(int k = jCol_[i]; k < jCol_[i + 1]; k++) {
2241        start[k] = nnz;
2242        length[k] = 0;
2243      }
2244    }
2245  }
2246  if(!needOrder) {
2247    {
2248      length[jCol_[nnz_jac_g -1] - 1]++;
2249      vals[nnz] = jValues_[nnz_jac_g - 1];
2250      rowLow[jRow_[nnz_jac_g - 1] - 1] += jValues_[nnz_jac_g - 1] * getColSolution()[jCol_ [nnz_jac_g - 1] - 1];
2251      rowUp[jRow_[nnz_jac_g - 1] -1] += jValues_[nnz_jac_g - 1] *getColSolution()[jCol_[nnz_jac_g - 1] -1];
2252      inds[nnz++] = jRow_[nnz_jac_g - 1] - 1;
2253    }
2254    for(int i = jCol_[nnz_jac_g -1] ; i < n ; i++) {
2255      start[i] = nnz;
2256      length[i] = 0;
2257    }
2258    start[n]=nnz;
2259    mat.assignMatrix(false, m, n, nnz, vals, inds, start, length);
2260  }
2261  else {
2262    std::cerr<<"jacobian matrix is not ordered properly"<<std::endl;
2263    throw -1;
2264  }
2265  delete [] g;
2266  int numcols=getNumCols();
2267  double *obj = new double[numcols];
2268  for(int i = 0 ; i < numcols ; i++)
2269    obj[i] = 0;
2270  mat.transpose();
2271#if 0
2272  std::cout<<"Mat is ordered by "<<
2273  <<mat.isColOrdered?"columns.":"rows."
2274  <<std::endl;
2275#endif
2276  si.loadProblem(mat, getColLower(), getColUpper(), obj, rowLow, rowUp);
2277  delete [] rowLow;
2278  delete [] rowUp;
2279  for(int i = 0 ; i < getNumCols() ; i++) {
2280    if(isInteger(i))
2281      si.setInteger(i);
2282  }
2283  if(getObj) {
2284    //add variable alpha
2285    //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
2286    CoinPackedVector a;
2287    si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
2288
2289    // Now get the objective cuts
2290    // get the gradient, pack it and add the cut
2291    tminlp_->eval_grad_f(n, getColSolution(), 1,obj);
2292    double ub;
2293    tminlp_->eval_f(n, getColSolution(), 1, ub);
2294    ub*=-1;
2295    double lb = -1e300;
2296    CoinPackedVector objCut;
2297    CoinPackedVector * v = &objCut;
2298    v->reserve(n);
2299    for(int i = 0; i<n ; i++) {
2300      if(cleanNnz(obj[i],colLower[i], colUpper[i],
2301          -getInfinity(), 0,
2302          getColSolution()[i],
2303          lb,
2304          ub, tiny_, veryTiny_)) {
2305        v->insert(i,obj[i]);
2306        lb += obj[i] * getColSolution()[i];
2307        ub += obj[i] * getColSolution()[i];
2308      }
2309    }
2310    v->insert(n,-1);
2311    si.addRow(objCut, lb, ub);
2312    delete [] obj;
2313  }
2314
2315// si.writeMps("init");
2316
2317  setWarmStartOptions();
2318  setColSolution(problem()->x_sol());
2319  setRowPrice(problem()->duals_sol());
2320
2321}
2322
2323
Note: See TracBrowser for help on using the repository browser.