source: branches/devel/Bonmin/src/Interfaces/BonOsiTMINLPInterface.cpp @ 101

Last change on this file since 101 was 101, checked in by pbonami, 13 years ago

Fixes a problem with random point generation with filter, and two small leaks

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