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

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

Add option to choose branching strategy,
include cmath in BonCurvatureEstimator?

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