source: trunk/Couenne/src/problem/CouenneProblem.cpp @ 78

Last change on this file since 78 was 78, checked in by pbelotti, 11 years ago

increase nOrigVars_ in CouenneProblem::addVariables (cleaner)

File size: 14.3 KB
Line 
1/*
2 * Name:    CouenneProblem.cpp
3 * Author:  Pietro Belotti
4 * Purpose: methods of the class CouenneProblem
5 *
6 * (C) Carnegie-Mellon University, 2006-09.
7 * This file is licensed under the Common Public License (CPL)
8 */
9
10#include <vector>
11
12#include "CoinHelperFunctions.hpp"
13#include "CoinTime.hpp"
14
15#include "CouenneTypes.hpp"
16
17#include "expression.hpp"
18#include "exprConst.hpp"
19#include "exprQuad.hpp"
20#include "exprClone.hpp"
21#include "exprIVar.hpp"
22#include "exprAux.hpp"
23#include "exprOpp.hpp"
24
25#include "CouenneProblem.hpp"
26#include "CouenneProblemElem.hpp"
27#include "depGraph.hpp"
28#include "lqelems.hpp"
29
30/// constructor
31CouenneProblem::CouenneProblem (const struct ASL *asl,
32                                Bonmin::BabSetupBase *base,
33                                JnlstPtr jnlst):
34  problemName_ (""),
35  auxSet_    (NULL), 
36  curnvars_  (-1),
37  nIntVars_  (0),
38  optimum_   (NULL),
39  bestObj_   (COIN_DBL_MAX),
40  quadIndex_ (NULL),
41  commuted_  (NULL),
42  numbering_ (NULL),
43  ndefined_  (0),
44  graph_     (NULL),
45  nOrigVars_ (0),
46  nOrigIntVars_ (0),
47  pcutoff_   (new GlobalCutOff (COIN_DBL_MAX)),
48  created_pcutoff_ (true),
49  doFBBT_    (true),
50  doOBBT_    (false),
51  doABT_     (true),
52  logObbtLev_(0),
53  logAbtLev_ (0),
54  jnlst_     (jnlst),
55  opt_window_ (COIN_DBL_MAX),
56  useQuadratic_ (false),
57  feas_tolerance_ (feas_tolerance_default),
58  integerRank_ (NULL),
59  maxCpuTime_  (COIN_DBL_MAX),
60  bonBase_     (base)
61{
62
63  double now = CoinCpuTime ();
64
65  if (asl) {
66#if COIN_HAS_ASL
67    // read problem from AMPL structure
68    readnl (asl);
69#else
70    jnlst_ -> Printf (Ipopt::J_ERROR, J_PROBLEM, "Couenne was compiled without ASL library. Cannot process ASL structure.\n");
71    throw -1;
72#endif
73
74    if ((now = (CoinCpuTime () - now)) > 10.)
75      jnlst_ -> Printf (Ipopt::J_WARNING, J_PROBLEM,
76                        "Couenne: reading time %.3fs\n", now);
77  }
78
79  // create expression set for binary search
80  auxSet_ = new std::set <exprAux *, compExpr>;
81
82  if (base) {
83    std::string s;
84    base -> options() -> GetStringValue ("use_quadratic", s, "couenne."); 
85    useQuadratic_ = (s == "yes");
86  }
87
88  if (base) {
89
90    std::string s;
91
92    base -> options() -> GetStringValue ("feasibility_bt",  s, "couenne."); doFBBT_ = (s == "yes");
93    base -> options() -> GetStringValue ("optimality_bt",   s, "couenne."); doOBBT_ = (s == "yes");
94    base -> options() -> GetStringValue ("aggressive_fbbt", s, "couenne."); doABT_  = (s == "yes");
95
96    base -> options() -> GetIntegerValue ("log_num_obbt_per_level", logObbtLev_, "couenne.");
97    base -> options() -> GetIntegerValue ("log_num_abt_per_level",  logAbtLev_,  "couenne.");
98
99    base -> options() -> GetNumericValue ("feas_tolerance",  feas_tolerance_, "couenne.");
100    base -> options() -> GetNumericValue ("opt_window",      opt_window_,     "couenne.");
101  }
102}
103
104
105/// preprocess problem in order to extract linear relaxations etc.
106void CouenneProblem::reformulate () {
107
108  double now = CoinCpuTime ();
109
110  if (domain_.current () == NULL) {
111
112    // create room for problem's variables and bounds, if no domain exists
113    CouNumber
114      *= (CouNumber *) malloc ((nVars()) * sizeof (CouNumber)),
115      *lb = (CouNumber *) malloc ((nVars()) * sizeof (CouNumber)),
116      *ub = (CouNumber *) malloc ((nVars()) * sizeof (CouNumber));
117
118    for (int i = nVars(); i--;) {
119      x  [i] =  0.;
120      lb [i] = -COUENNE_INFINITY;
121      ub [i] =  COUENNE_INFINITY;
122    }
123
124    domain_.push (nVars (), x, lb, ub);
125  }
126
127  // link initial variables to problem's domain
128  for (std::vector <exprVar *>::iterator i = variables_.begin ();
129       i != variables_.end (); ++i)
130    (*i) -> linkDomain (&domain_);
131
132  if (jnlst_ -> ProduceOutput(Ipopt::J_SUMMARY, J_PROBLEM))
133    print (std::cout);
134
135  // save -- for statistic purposes -- number of original
136  // constraints. Some of them will be deleted as definition of
137  // auxiliary variables.
138  nOrigCons_    = constraints_. size ();
139  nOrigIntVars_ = nIntVars ();
140
141  jnlst_->Printf (Ipopt::J_SUMMARY, J_PROBLEM,
142                  "Problem size before reformulation: %d variables (%d integer), %d constraints.\n",
143                  nOrigVars_, nOrigIntVars_, nOrigCons_);
144
145  // reformulation
146  if (!standardize ()) { // problem is infeasible if standardize returns false
147
148    jnlst_->Printf(Ipopt::J_WARNING, J_PROBLEM,
149                   "Couenne: problem infeasible after reformulation\n");
150    // fake infeasible bounds for Couenne to bail out
151    for (int i = nVars (); i--;)
152      Ub (i) = - (Lb (i) = 1.);
153  }
154
155  // clear all spurious variables pointers not referring to the variables_ vector
156  realign ();
157
158  // give a value to all auxiliary variables. Do it now to be able to
159  // recognize complementarity constraints in fillDependence()
160  initAuxs ();
161
162  // fill dependence_ structure
163  fillDependence (bonBase_);
164
165  // quadratic handling
166  fillQuadIndices ();
167
168  if ((now = (CoinCpuTime () - now)) > 10.)
169    jnlst_->Printf(Ipopt::J_WARNING, J_PROBLEM,
170    "Couenne: reformulation time %.3fs\n", now);
171
172  // give a value to all auxiliary variables
173  initAuxs ();
174
175  int nActualVars = nIntVars_ = 0;
176
177  // check how many integer variables we have now (including aux)
178  for (int i=0; i<nVars(); i++)
179    if (variables_ [i] -> Multiplicity () > 0) {
180
181      nActualVars++;
182      if (variables_ [i] -> isDefinedInteger ())
183        nIntVars_++;
184    }
185
186  jnlst_->Printf(Ipopt::J_SUMMARY, J_PROBLEM,
187                  "Problem size after  reformulation: %d variables (%d integer), %d constraints.\n",
188                  nActualVars, nIntVars_, nCons());
189
190  // check if optimal solution is available (for debug purposes)
191  readOptimum ();
192
193  if (bonBase_) {
194
195    CouNumber
196      art_cutoff =  COIN_DBL_MAX,
197      art_lower  = -COIN_DBL_MAX;
198
199    bonBase_ -> options() -> GetNumericValue ("art_cutoff", art_cutoff, "couenne.");
200    bonBase_ -> options() -> GetNumericValue ("art_lower",  art_lower,  "couenne.");
201
202    if (art_cutoff <  1.e50) setCutOff (art_cutoff);
203    if (art_lower  > -1.e50) {
204      int indobj = objectives_ [0] -> Body () -> Index ();
205      if (indobj >= 0)
206        domain_.lb (indobj) = art_lower;
207    }
208  }
209
210  if (jnlst_->ProduceOutput(Ipopt::J_DETAILED, J_PROBLEM)) {
211    // We should route that also through the journalist
212    print (std::cout);
213  }
214
215  //writeAMPL ("extended-aw.mod", true);
216  //writeAMPL ("original.mod", false);
217}
218
219
220/// copy constructor
221
222CouenneProblem::CouenneProblem (const CouenneProblem &p):
223  problemName_  (p.problemName_),
224  domain_       (p.domain_),
225  curnvars_     (-1),
226  nIntVars_     (p.nIntVars_),
227  optimum_      (NULL),
228  bestObj_      (p.bestObj_),
229  commuted_     (NULL),
230  numbering_    (NULL),
231  ndefined_     (p.ndefined_),
232  graph_        (NULL),
233  nOrigVars_    (p.nOrigVars_),
234  nOrigCons_    (p.nOrigCons_),
235  nOrigIntVars_ (p.nOrigIntVars_),
236  pcutoff_      (p.pcutoff_),
237  created_pcutoff_ (false),
238  doFBBT_       (p. doFBBT_),
239  doOBBT_       (p. doOBBT_),
240  doABT_        (p. doABT_),
241  logObbtLev_   (p. logObbtLev_),
242  logAbtLev_    (p. logAbtLev_),
243  jnlst_        (p.jnlst_),
244  opt_window_   (p.opt_window_),    // needed only in standardize (), unnecessary to update it
245  useQuadratic_ (p.useQuadratic_),  // ditto
246  feas_tolerance_ (p.feas_tolerance_),
247  dependence_   (p.dependence_),
248  objects_      (p.objects_),
249  integerRank_  (NULL),
250  numberInRank_ (p.numberInRank_),
251  maxCpuTime_   (p.maxCpuTime_),
252  bonBase_      (p.bonBase_) {
253
254  for (int i=0; i < p.nVars (); i++)
255    variables_ . push_back (NULL);
256
257  for (int i=0; i < p.nVars (); i++) {
258    int ind = p.numbering_ [i];
259    variables_ [ind] = p.Var (ind) -> clone (&domain_);
260  }
261
262  if (p.numbering_)
263    numbering_ = CoinCopyOfArray (p.numbering_, nVars ());
264
265  // clone objectives and constraints (there's a leak around here)
266  for (int i=0; i < p.nObjs (); i++) objectives_  . push_back (p.Obj (i) -> clone (&domain_));
267  for (int i=0; i < p.nCons (); i++) constraints_ . push_back (p.Con (i) -> clone (&domain_));
268
269  if (p.optimum_) 
270    optimum_ = CoinCopyOfArray (p.optimum_, nVars ());
271   
272  // clear all spurious variables pointers not referring to the variables_ vector
273  realign ();
274
275  // copy integer rank (used in getIntegerCandidate)
276  if (p.integerRank_) {
277    integerRank_ = new int [nVars ()];
278    CoinCopyN (p.integerRank_, nVars (), integerRank_);
279  }
280}
281
282
283/// Destructor
284
285CouenneProblem::~CouenneProblem () {
286
287  // delete optimal solution (if any)
288  if (optimum_)
289    free (optimum_);
290
291  // delete objectives
292  for (std::vector <CouenneObjective *>::iterator i  = objectives_ . begin ();
293       i != objectives_ . end (); ++i)
294    delete (*i);
295
296  // delete constraints
297  for (std::vector <CouenneConstraint *>::iterator i = constraints_ . begin ();
298       i != constraints_ . end (); ++i)
299    delete (*i);
300
301  // delete variables
302  //for (std::vector <exprVar *>::iterator i = variables_ . begin ();
303  //i != variables_ . end (); ++i)
304  //delete (*i);
305
306  for (int i=nVars (); i--;) { // delete in inverse order
307    int ind = numbering_ [i];
308    delete variables_ [ind];
309  }
310
311  // delete extra structures
312  if (graph_)     delete    graph_;
313  if (commuted_)  delete [] commuted_;
314  if (numbering_) delete [] numbering_;
315
316  if (created_pcutoff_) delete pcutoff_;
317
318  if (integerRank_) delete [] integerRank_;
319}
320
321
322/// methods to add objective function.
323
324void CouenneProblem::addObjective (expression *newobj, const std::string &sense = "min") {
325  objectives_ . push_back
326    (new CouenneObjective ((sense == "min") ? newobj : new exprOpp (new exprClone (newobj))));//, MINIMIZE));
327}
328
329
330/// methods to add nonlinear constraints:
331
332/// equality constraint
333void CouenneProblem::addEQConstraint (expression *body, expression *rhs = NULL) {
334
335  if (!rhs) rhs = new exprConst (0.);
336  constraints_ . push_back (new CouenneConstraint (body, rhs, new exprClone (rhs)));
337}
338
339/// "greater than" constraint
340void CouenneProblem::addGEConstraint (expression *body, expression *rhs = NULL) {
341  if (!rhs) rhs = new exprConst (0.);
342  constraints_ . push_back (new CouenneConstraint
343                            (body, rhs, new exprConst (COUENNE_INFINITY)));
344}
345
346/// "smaller than" constraint
347void CouenneProblem::addLEConstraint (expression *body, expression *rhs = NULL) {
348  if (!rhs) rhs = new exprConst (0.);
349  constraints_ . push_back (new CouenneConstraint
350                            (body, new exprConst (-COUENNE_INFINITY), rhs));
351}
352
353/// range constraint
354void CouenneProblem::addRNGConstraint (expression *body, expression *lb=NULL, expression *ub=NULL) {
355  if (!lb) lb = new exprConst (0.);
356  if (!ub) ub = new exprConst (0.);
357  constraints_ . push_back (new CouenneConstraint (body, lb, ub));
358}
359
360
361/// add variable to the problem -- check whether it is integer (isDiscrete)
362
363expression *CouenneProblem::addVariable (bool isDiscrete, Domain *d) {
364
365  exprVar *var = (isDiscrete) ? 
366    (new exprIVar (variables_ . size (), d)) :
367    (new exprVar  (variables_ . size (), d));
368
369  variables_ . push_back (var);
370
371  if (isDiscrete) 
372    nIntVars_++;
373
374  nOrigVars_++;
375
376  return var;
377}
378
379
380/// add auxiliary variable and associate it with pointer to expression
381/// given as argument
382exprAux *CouenneProblem::addAuxiliary (expression *symbolic) {
383
384  // check if image is already in the expression database auxSet_
385  std::set <exprAux *, compExpr>::iterator i;
386
387  int var_ind = variables_ . size ();
388  domain_. current () -> resize (var_ind + 1);
389
390  symbolic -> getBounds (domain_. lb (var_ind), 
391                         domain_. ub (var_ind));
392
393  // create new aux associated with that expression
394  exprAux *w = new exprAux (symbolic,
395                            var_ind,
396                            1 + symbolic -> rank (), 
397                            exprAux::Unset, 
398                            &domain_);
399  //symbolic -> isInteger () ? exprAux::Integer : exprAux::Continuous);
400
401  //  w -> linkDomain (&domain_);
402
403  // seek expression in the set
404  if ((i = auxSet_ -> find (w)) == auxSet_ -> end ()) {
405
406    // no such expression found in the set, create entry therein
407    variables_ . push_back (w);
408    auxSet_ -> insert (w); // insert into repetition checking structure
409    graph_  -> insert (w); // insert into acyclic structure
410
411  } else {  // otherwise, just return the entry's pointer
412
413    delete w;
414    w = *i;
415    (*i) -> increaseMult ();
416  }
417
418  return w;
419}
420
421
422/// translates pair (indices, coefficients) into vector with pointers to variables
423void CouenneProblem::indcoe2vector (int *indexL, 
424                                    CouNumber *coeff,
425                                    std::vector <std::pair <exprVar *, CouNumber> > &lcoeff) {
426  // TODO: sort
427
428  for (int i=0; indexL [i] >= 0; i++)
429    lcoeff.push_back (std::pair <exprVar *, CouNumber> (Var (indexL [i]), coeff [i]));
430}
431
432
433/// translates triplet (indicesI, indicesJ, coefficients) into vector with pointers to variables
434void CouenneProblem::indcoe2vector (int *indexI,
435                                    int *indexJ,
436                                    CouNumber *coeff,
437                                    std::vector <quadElem> &qcoeff) {
438  // TODO: sort
439
440  for (int i=0; indexI [i] >= 0; i++)
441    qcoeff.push_back (quadElem (Var (indexI [i]), Var (indexJ [i]), coeff [i]));
442}
443
444
445/// fill in the integerRank_ array
446void CouenneProblem::fillIntegerRank () const {
447
448  if (integerRank_)
449    return;
450
451  int nvars = nVars ();
452
453  integerRank_ = new int [nvars];
454
455  // 0: fractional
456  // 1: integer
457  // k: integer,    depending on at least one integer with associated value k-1, or
458  // k: fractional, depending on at least one integer with associated value k
459
460  for (int ii = 0; ii < nvars; ii++) {
461
462    int index = numbering_ [ii];
463
464    if (Var (index) -> Multiplicity () <= 0) {
465      integerRank_ [index] = 0;
466      continue;
467    }
468
469    bool isInt = Var (index) -> isDefinedInteger ();
470
471    integerRank_ [index] = (isInt) ? 1 : 0;
472
473    if (Var (index) -> Type () == AUX) {
474
475      std::set <int> deplist;
476
477      if (Var (index) -> Image () -> DepList (deplist, STOP_AT_AUX) != 0) // depends on something
478        for (std::set <int>::iterator i = deplist.begin (); i != deplist.end (); ++i) {
479
480          int token = integerRank_ [*i];
481          if (isInt) token++;
482
483          if (token > integerRank_ [index]) // there's a free integer below us
484            integerRank_ [index] = token;
485        }
486    }
487  }
488
489  jnlst_->Printf (Ipopt::J_VECTOR, J_PROBLEM, "Free (original) integers\n");
490  for (int i=0; i<nOrigVars_; i++)
491    jnlst_->Printf (Ipopt::J_VECTOR, J_PROBLEM, "%d: %d\n", i, integerRank_ [i]);
492
493  // fill in numberInRank_
494  for (int i=0; i<nOrigVars_; i++)
495    if ((variables_ [i] -> isDefinedInteger ()) &&
496        (variables_ [i] -> Multiplicity () > 0)) {
497
498    int rank = integerRank_ [i];
499
500    if (numberInRank_.size () <= (unsigned int) rank)
501      for (int j=numberInRank_.size (); j <= rank; j++)
502        numberInRank_ .push_back (0);
503
504    numberInRank_ [rank] ++;
505  }
506
507  jnlst_->Printf (Ipopt::J_VECTOR, J_PROBLEM, "numInteger [neglect non-originals]\n");
508  for (unsigned int i=0; i<numberInRank_.size(); i++)
509    jnlst_->Printf (Ipopt::J_VECTOR, J_PROBLEM, "%d: %d\n", i, numberInRank_ [i]);
510}
Note: See TracBrowser for help on using the repository browser.