source: branches/devel/Bonmin/src/Interfaces/BonTMINLP2TNLP.cpp @ 104

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

Fix a bug with zero dimensional nlps

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 16.2 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// Carl D. Laird, Carnegie Mellon University,
7// Andreas Waechter, International Business Machines Corporation
8// Pierre Bonami, Carnegie Mellon University,
9//
10// Date : 12/01/2004
11
12
13#include "BonTMINLP2TNLP.hpp"
14#include "IpBlas.hpp"
15#include "IpAlgTypes.hpp"
16#include <string>
17#include <fstream>
18#include <sstream>
19namespace Bonmin
20{
21  TMINLP2TNLP::TMINLP2TNLP(const SmartPtr<TMINLP> tminlp
22#ifdef WARM_STARTER
23       ,
24      const OptionsList& options
25#endif
26       )
27      :
28      tminlp_(tminlp),
29      n_(0),
30      m_(0),
31      var_types_(NULL),
32      x_l_(NULL),
33      x_u_(NULL),
34      orig_x_l_(NULL),
35      orig_x_u_(NULL),
36      g_l_(NULL),
37      g_u_(NULL),
38      x_init_(NULL),
39      x_init_user_(NULL),
40      x_sol_(NULL),
41      g_sol_(NULL),
42      duals_sol_(NULL),
43//      return_status_(NOT_SOLVED),
44      obj_value_(1e100),
45      curr_warm_starter_(NULL),
46      need_new_warm_starter_(true)
47  {
48    // read the nlp size and bounds information from
49    // the TMINLP and keep an internal copy. This way the
50    // caller can modify the bounds that are sent to Ipopt;
51    DBG_ASSERT(IsValid(tminlp_));
52
53    Index nnz_jac_g;
54    tminlp_->get_nlp_info(n_, m_, nnz_jac_g, nnz_h_lag_, index_style_);
55
56    // Allocate space for the variable types vector
57    var_types_ = new TMINLP::VariableType[n_];
58
59    // retrieve the variable types
60    tminlp_->get_variables_types(n_, var_types_);
61
62    // Allocate space for the internal copy of the variable bounds
63    x_l_ = new Number[n_];
64    x_u_ = new Number[n_];
65    orig_x_l_ = new Number[n_];
66    orig_x_u_ = new Number[n_];
67
68    g_l_ = new Number[m_];
69    g_u_ = new Number[m_];
70
71    // retrieve the variable bounds
72    tminlp_->get_bounds_info(n_, x_l_, x_u_, m_, g_l_, g_u_);
73    IpBlasDcopy(n_, x_l_, 1, orig_x_l_, 1);
74    IpBlasDcopy(n_, x_u_, 1, orig_x_u_, 1);
75
76    //    // Check that the bounds make sense compared with the variable type
77    //    for (int i=0; i<n_; i++) {
78    //      throw_exception_on_bad_variable_bound(i);
79    //    }
80
81    // Allocate space for the initial point
82    x_init_ = new Number[3*n_ + m_];
83    tminlp_->get_starting_point(n_, true, x_init_, false, NULL, NULL,
84        m_, false, NULL);
85    x_init_user_ = new Number[n_];
86    IpBlasDcopy(n_, x_init_, 1, x_init_user_, 1);
87    duals_sol_ = NULL;
88    duals_init_ = NULL;
89
90#ifdef WARM_STARTER
91    // Get values for parameters
92    options.GetNumericValue("nlp_lower_bound_inf", nlp_lower_bound_inf_, "");
93    options.GetNumericValue("nlp_upper_bound_inf", nlp_upper_bound_inf_, "");
94    options.GetBoolValue("warm_start_entire_iterate",
95        warm_start_entire_iterate_, "");
96#endif
97  }
98
99  TMINLP2TNLP::~TMINLP2TNLP()
100  {
101    delete [] var_types_;
102    var_types_ = NULL;
103    delete [] x_l_;
104    x_l_ = NULL;
105    delete [] x_u_;
106    x_u_ = NULL;
107    delete [] orig_x_l_;
108    orig_x_l_ = NULL;
109    delete [] orig_x_u_;
110    orig_x_u_ = NULL;
111    delete [] g_l_;
112    g_l_ = NULL;
113    delete [] g_u_;
114    g_u_ = NULL;
115    delete [] x_init_;
116    x_init_ = NULL;
117    delete [] x_init_user_;
118    x_init_user_ = NULL;
119    delete [] x_sol_;
120    x_sol_ = NULL;
121    delete [] g_sol_;
122    g_sol_ = NULL;
123    delete [] duals_sol_;
124    duals_sol_ = NULL;
125  }
126
127  /** Copies the modification made to TNLP by the user (modifications
128      such as changing bound changing starting point,...).
129      this and other should be two instances of the same problem
130      I am trying to mimic a copy construction for Cbc
131      use with great care not safe.
132  */
133  void
134  TMINLP2TNLP::copyUserModification(TMINLP2TNLP& other)
135  {
136    DBG_ASSERT(x_l_);
137    DBG_ASSERT(x_u_);
138    DBG_ASSERT(other.x_l_);
139    DBG_ASSERT(other.x_u_);
140    DBG_ASSERT(n_ == other.n_);
141    DBG_ASSERT(m_ == other.m_);
142
143    IpBlasDcopy(n_, other.x_l_, 1, x_l_, 1);
144    IpBlasDcopy(n_, other.x_u_, 1, x_u_, 1);
145
146    if(other.duals_init_) {
147      duals_init_ = &x_init_[n_];
148      IpBlasDcopy(3*n_+ m_, other.x_init_, 1, x_init_, 1);
149    }
150    else
151      IpBlasDcopy(n_, other.x_init_, 1, x_init_, 1);
152    return_status_ = other.return_status_;
153
154    if(other.x_sol_ !=NULL) {
155      //DBG_ASSERT(return_status_ != NOT_SOLVED);
156      Set_x_sol(n_,other.x_sol_);
157    }
158
159    if(other.g_sol_!=NULL) {
160//      DBG_ASSERT(return_status_ != NOT_SOLVED);
161      g_sol_ = new Number [m_];
162      IpBlasDcopy(m_, other.g_sol_, 1, g_sol_, 1);
163    }
164
165    if(other.duals_sol_!=NULL) {
166//      DBG_ASSERT(return_status_ != NOT_SOLVED);
167      duals_sol_ = new Number[m_ + 2*n_];
168      IpBlasDcopy(2*n_+ m_, other.duals_sol_, 1, duals_sol_, 1);
169    }
170
171    obj_value_ = other.obj_value_;
172
173    curr_warm_starter_ = other.curr_warm_starter_;
174
175    nlp_lower_bound_inf_ = other.nlp_lower_bound_inf_;
176
177    nlp_upper_bound_inf_ = other.nlp_upper_bound_inf_;
178
179    need_new_warm_starter_ = other.need_new_warm_starter_;
180  }
181
182  void TMINLP2TNLP::SetVariableBounds(Index var_no, Number x_l, Number x_u)
183  {
184    DBG_ASSERT(var_no >= 0 && var_no < n_);
185    x_l_[var_no] = x_l;
186    x_u_[var_no] = x_u;
187    //    throw_exception_on_bad_variable_bound(var_no);
188  }
189
190  void TMINLP2TNLP::SetVariableLowerBound(Index var_no, Number x_l)
191  {
192    DBG_ASSERT(var_no >= 0 && var_no < n_);
193    x_l_[var_no] = x_l;
194    //    throw_exception_on_bad_variable_bound(var_no);
195  }
196
197  void TMINLP2TNLP::SetVariableUpperBound(Index var_no, Number x_u)
198  {
199    DBG_ASSERT(var_no >= 0 && var_no < n_);
200    x_u_[var_no] = x_u;
201    //    throw_exception_on_bad_variable_bound(var_no);
202  }
203
204  void TMINLP2TNLP::SetStartingPoint(Index n,const Number* x_init)
205  {
206    DBG_ASSERT(n == n_);
207    IpBlasDcopy(n_, x_init, 1, x_init_, 1);
208  }
209
210  void TMINLP2TNLP::resetStartingPoint()
211  {
212    curr_warm_starter_ = NULL;
213    IpBlasDcopy(n_, x_init_user_, 1, x_init_, 1);
214  }
215
216  void TMINLP2TNLP::setxInit(Index ind, const Number val)
217  {
218    x_init_[ind] = val;
219  }
220
221  void TMINLP2TNLP::setxInit(Index n,const Number* x_init)
222  {
223    DBG_ASSERT(n == n_);
224    IpBlasDcopy(n_, x_init, 1, x_init_, 1);
225  }
226
227  void TMINLP2TNLP::setDualInit(Index ind, const Number val)
228  {
229    if(!duals_init_)
230      duals_init_ = &x_init_[n_];
231    duals_init_[ind] = val;
232  }
233
234  void TMINLP2TNLP::setDualsInit(Index m, const Number* duals_init)
235  {
236    DBG_ASSERT(m == m_ + 2*n_ );
237
238    if(!duals_init_)
239      duals_init_ = x_init_ + n_;
240
241    IpBlasDcopy(m, duals_init, 1, duals_init_, 1);
242
243  }
244
245  /** Set the contiuous solution */
246  void TMINLP2TNLP::Set_x_sol(Index n, const Number* x_sol)
247  {
248    DBG_ASSERT(n == n_);
249    if (!x_sol_) {
250      x_sol_ = new Number[n];
251    }
252    IpBlasDcopy(n_, x_sol, 1, x_sol_, 1);
253  }
254
255  /** Change the type of the variable */
256  void TMINLP2TNLP::SetVariableType(Index n, TMINLP::VariableType type)
257  {
258    DBG_ASSERT(n >= 0 && n < n_);
259    var_types_[n] = type;
260  }
261
262  bool TMINLP2TNLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
263      Index& nnz_h_lag, TNLP::IndexStyleEnum& index_style)
264  {
265    bool return_code = tminlp_->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
266    m += tminlp_->nLinearCuts_;
267    nnz_jac_g += tminlp_->linearCutsNnz_;
268    return return_code;
269  }
270
271  bool TMINLP2TNLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
272      Index m, Number* g_l, Number* g_u)
273  {
274    DBG_ASSERT(n==n_);
275    DBG_ASSERT(m==m_);
276    IpBlasDcopy(n_, x_l_, 1, x_l, 1);
277    IpBlasDcopy(n_, x_u_, 1, x_u, 1);
278    IpBlasDcopy(m_, g_l_, 1, g_l, 1);
279    IpBlasDcopy(m_, g_u_, 1, g_u, 1);
280    IpBlasDcopy(tminlp_->nLinearCuts_, tminlp_->lower_, 1, &g_l[m_],1); 
281    IpBlasDcopy(tminlp_->nLinearCuts_, tminlp_->upper_, 1, &g_u[m_],1); 
282    return true;
283  }
284
285  bool TMINLP2TNLP::get_starting_point(Index n, bool init_x, Number* x,
286      bool init_z, Number* z_L, Number* z_U,
287      Index m, bool init_lambda,
288      Number* lambda)
289  {
290    assert(m==m_ + tminlp_->nLinearCuts_);
291    assert(n==n_);
292    if (init_x == true) {
293      if(x_init_==NULL)
294        return false;
295      IpBlasDcopy(n, x_init_, 1, x, 1);
296    }
297    if (init_z == true) {
298      if(duals_init_ == NULL)
299        return false;
300      IpBlasDcopy(n, &duals_init_[m], 1, z_L, 1);
301      IpBlasDcopy(n, &duals_init_[m + n], 1, z_U, 1);
302
303    }
304    if(init_lambda == true) {
305      if(duals_init_ == NULL)
306        return false;
307      IpBlasDcopy(m_, duals_init_, 1, lambda, 1);
308      for(int i = 0 ; i < tminlp_->nLinearCuts_; i++)
309      {
310        lambda [i + m_ ] = 0.;
311      }
312    }
313
314    need_new_warm_starter_ = true;
315    return true;
316  }
317
318  bool TMINLP2TNLP::get_warm_start_iterate(IteratesVector& warm_start_iterate)
319  {
320    if (IsNull(curr_warm_starter_)) {
321      return false;
322    }
323
324    bool retval = curr_warm_starter_->WarmStartIterate(n_, x_l_, x_u_,
325        warm_start_iterate);
326
327    need_new_warm_starter_ = true;
328    return retval;
329  }
330
331  bool TMINLP2TNLP::eval_f(Index n, const Number* x, bool new_x,
332      Number& obj_value)
333  {
334    return tminlp_->eval_f(n, x, new_x, obj_value);
335  }
336
337  bool TMINLP2TNLP::eval_grad_f(Index n, const Number* x, bool new_x,
338      Number* grad_f)
339  {
340    return tminlp_->eval_grad_f(n, x, new_x, grad_f);
341  }
342
343  bool TMINLP2TNLP::eval_g(Index n, const Number* x, bool new_x,
344      Index m, Number* g)
345  {
346    int return_code = tminlp_->eval_g(n, x, new_x, m_, g);
347    // Add the linear cuts
348    int nnz = 0;
349    for(int i = 0 ; i < tminlp_->nLinearCuts_ ; i++)
350    {
351      int iplusm = m_ + i;
352      g[iplusm] = 0;
353      while(tminlp_->iRow_[nnz]==i) { g[iplusm] += tminlp_->elems_[nnz] * x[tminlp_->jCol_[nnz]]; nnz++;}
354    }
355    return return_code;
356  }
357
358  bool TMINLP2TNLP::eval_jac_g(Index n, const Number* x, bool new_x,
359      Index m, Index nele_jac, Index* iRow,
360      Index *jCol, Number* values)
361  {
362    int return_code = tminlp_->eval_jac_g(n, x, new_x, m_ , nele_jac - tminlp_->nLinearCuts_, 
363                                          iRow, jCol, values);
364    bool isFortran = index_style_ == TNLP::FORTRAN_STYLE;
365    if(iRow != NULL)
366    {
367      DBG_ASSERT(jCol != NULL);
368      DBG_ASSERT(values == NULL);
369     
370      int nnz = nele_jac - tminlp_->linearCutsNnz_ ;
371      for(int i = 0; i < tminlp_->linearCutsNnz_ ; i++ , nnz++)
372      {
373         iRow[nnz] = tminlp_->iRow_[i] + m_ + isFortran; 
374         jCol[nnz] = tminlp_->jCol_[i] + isFortran;
375      }
376    }
377    else
378    {
379      DBG_ASSERT(jCol == NULL);
380      DBG_ASSERT(values != NULL);
381      IpBlasDcopy(tminlp_->linearCutsNnz_ , tminlp_->elems_, 1, &values[nele_jac - tminlp_->linearCutsNnz_],1);
382    }
383  return return_code;
384  }
385
386  bool TMINLP2TNLP::eval_h(Index n, const Number* x, bool new_x,
387      Number obj_factor, Index m, const Number* lambda,
388      bool new_lambda, Index nele_hess,
389      Index* iRow, Index* jCol, Number* values)
390  {
391    return tminlp_->eval_h(n, x, new_x, obj_factor, m, lambda,
392        new_lambda, nele_hess,
393        iRow, jCol, values);
394  }
395
396  void TMINLP2TNLP::finalize_solution(SolverReturn status,
397      Index n, const Number* x, const Number* z_L, const Number* z_U,
398      Index m, const Number* g, const Number* lambda,
399      Number obj_value)
400  {
401    assert(n==n_);
402    assert(m==m_ + tminlp_->nLinearCuts_);
403    if (!x_sol_) {
404      x_sol_ = new Number[n];
405    }
406    IpBlasDcopy(n, x, 1, x_sol_, 1);
407
408    if(!g_sol_) {
409      g_sol_ = new Number [m];
410    }
411    IpBlasDcopy(m, g, 1, g_sol_, 1);
412    if (!duals_sol_) {
413      duals_sol_ = new Number[m + 2*n];
414    }
415    IpBlasDcopy(m, lambda, 1, duals_sol_, 1);
416
417    IpBlasDcopy(n, z_L, 1 , duals_sol_ + m, 1);
418    IpBlasDcopy(n, z_U, 1 , duals_sol_ + m + n, 1);
419
420    return_status_ = status;
421    obj_value_ = obj_value;
422
423    if (IsValid(curr_warm_starter_)) {
424      curr_warm_starter_->Finalize();
425    }
426  }
427
428  void TMINLP2TNLP::throw_exception_on_bad_variable_bound(Index i)
429  {
430    DBG_ASSERT(i >= 0 && i < n_);
431
432    if (var_types_[i] == TMINLP::BINARY) {
433      ASSERT_EXCEPTION(x_l_[i] == 0 && x_u_[i] == 1,
434          TMINLP_INVALID_VARIABLE_BOUNDS,
435          "Invalid variable bounds in TMINLP. All binaries must have 0,1 for variable bounds."
436                      );
437    }
438    else if (var_types_[i] == TMINLP::INTEGER) {
439      if(fabs(x_l_[i])<INT_MAX && fabs(x_u_[i]) < INT_MAX)//round to closest valid integer
440      {
441        int x_l = (int)ceil(x_l_[i]);
442        int x_u = (int)floor(x_u_[i]);
443        std::cerr<<"Inconsistent bounds on an integer"<<std::endl;
444        ASSERT_EXCEPTION(x_l_[i] == (Number)x_l && x_u_[i] == (Number)x_u,
445            TMINLP_INVALID_VARIABLE_BOUNDS,
446            "Invalid variable bounds in TMINLP. All integer variables must have integer bounds."
447                        );
448      }
449    }
450  }
451
452
453  bool TMINLP2TNLP::intermediate_callback(AlgorithmMode mode,
454      Index iter, Number obj_value,
455      Number inf_pr, Number inf_du,
456      Number mu, Number d_norm,
457      Number regularization_size,
458      Number alpha_du, Number alpha_pr,
459      Index ls_trials,
460      const IpoptData* ip_data,
461      IpoptCalculatedQuantities* ip_cq)
462  {
463#if WARM_STARTER
464    // If we don't have this swtiched on, we assume that also the
465    // "warm_start" option for bonmin is set not to refer to the
466    // interior warm start object
467    if (!warm_start_entire_iterate_) {
468      return true;
469    }
470    if (need_new_warm_starter_) {
471      // Create a new object for later warm start information
472      curr_warm_starter_ = new IpoptInteriorWarmStarter(n_, x_l_, x_u_,
473          nlp_lower_bound_inf_,
474          nlp_upper_bound_inf_,
475          warm_start_entire_iterate_);
476      need_new_warm_starter_ = false;
477    }
478
479    return curr_warm_starter_->UpdateStoredIterates(mode, *ip_data, *ip_cq);
480#else
481    return true;
482#endif
483  }
484
485
486  /** Procedure to ouptut relevant informations to reproduce a sub-problem.
487  Compare the current problem to the problem to solve
488  and writes files with bounds which have changed and current starting point.
489
490  */
491  void
492  TMINLP2TNLP::outputDiffs(const std::string& probName, const std::string * varNames)
493  {
494    const int &numcols = n_;
495    const int &numrows = m_;
496
497    const double * currentLower = x_l();
498    const double * currentUpper = x_u();
499
500    const double * originalLower = orig_x_l();
501    const double * originalUpper = orig_x_u();
502    CoinRelFltEq eq;
503    std::string fBoundsName = probName;
504    std::ostringstream os;
505    fBoundsName+=".bounds";
506    std::string fModName = probName;
507    fModName+= ".mod";
508    std::ofstream fBounds;
509    std::ofstream fMod;
510    bool hasVarNames = 0;
511
512    if(varNames!=NULL )
513      hasVarNames=1;
514    if(hasVarNames)
515      fMod.open(fModName.c_str());
516    fBounds.open(fBoundsName.c_str());
517
518    for(int i = 0 ; i < numcols ; i++) {
519      if(!eq(currentLower[i],originalLower[i])) {
520        if(hasVarNames)
521          fMod<<"bounds"<<i<<": "
522          <<varNames[i]<<" >= "
523          <<currentLower[i]<<";\n";
524
525
526        fBounds<<"LO"<<"\t"<<i<<"\t"<<currentLower[i]<<std::endl;
527      }
528      if(!eq(currentUpper[i],originalUpper[i])) {
529        if(hasVarNames)
530          fMod<<"bounds"<<i<<": "
531          <<varNames[i]<<" <= "
532          <<currentUpper[i]<<";\n";
533
534        fBounds<<"UP"<<"\t"<<i<<"\t"<<currentUpper[i]<<std::endl;
535      }
536    }
537
538    //write a file with starting point
539    std::string fStartPointName=probName;
540    fStartPointName+=".start";
541
542    std::ofstream fStartPoint(fStartPointName.c_str());
543    const double * primals = x_init();
544    const double * duals = duals_init();
545    fStartPoint.precision(17);
546    fStartPoint<<numcols<<"\t"<<2*numcols+numrows<<std::endl;
547    for(int i = 0 ; i < numcols ; i++)
548      fStartPoint<<primals[i]<<std::endl;
549    int end = 2*numcols + numrows;
550    if(duals) {
551      for(int i = 0 ; i < end; i++)
552        fStartPoint<<duals[i]<<std::endl;
553    }
554
555  }
556
557  /** force solution to be fractionnal.*/
558  void
559  TMINLP2TNLP::force_fractionnal_sol()
560  {
561    for(int i=0 ; i < n_ ; i++) {
562      if( ( var_types_[i] == TMINLP::INTEGER ||
563          var_types_[i] == TMINLP::BINARY )&&
564          x_l_[i] < x_u_[i] + 0.5)//not fixed
565      {
566        x_sol_[i] = ceil(x_l_[i]) + 0.5;//make it integer infeasible
567      }
568    }
569  }
570
571  bool 
572  TMINLP2TNLP::get_scaling_parameters(Number& obj_scaling,
573                                        bool& use_x_scaling, Index n,
574                                        Number* x_scaling,
575                                        bool& use_g_scaling, Index m,
576                                        Number* g_scaling)
577  {
578    return tminlp_->get_scaling_parameters(obj_scaling, use_x_scaling, n,
579                                  x_scaling,
580                                  use_g_scaling, m, g_scaling);
581  }
582                                 
583
584
585
586}
587// namespace Ipopt
588
Note: See TracBrowser for help on using the repository browser.