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

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

Resolve small problems in unitTest

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