source: trunk/Couenne/src/problem/problem.cpp @ 61

Last change on this file since 61 was 61, checked in by pbelotti, 12 years ago

fixed some method for quadratic functions. Deleted useless problem print.

File size: 11.4 KB
Line 
1/*
2 * Name:    problem.cpp
3 * Author:  Pietro Belotti
4 * Purpose: methods of the class CouenneProblem
5 *
6 * (C) Carnegie-Mellon University, 2006-08.
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 "exprGroup.hpp"
20#include "exprClone.hpp"
21#include "exprAux.hpp"
22#include "lqelems.hpp"
23#include "CouenneProblem.hpp"
24#include "CouenneProblemElem.hpp"
25#include "lqelems.hpp"
26
27const CouNumber SafeCutoff = COUENNE_EPS;
28
29/// initialize auxiliary variables from original variables in the
30/// nonlinear problem
31
32void CouenneProblem::initAuxs () const {
33
34  domain_.current () -> resize (nVars ());
35
36  // initially, auxiliary variables are unbounded, their bounds only
37  // depending on their function
38
39  int nvars = nVars ();
40
41  for (int i=0; i < nvars; i++) {
42
43    int indvar = variables_ [i] -> Index ();
44
45    if ((variables_ [i] -> Type () == AUX) &&                   // this is an auxiliary
46        (indvar >= nOrigVars_) || // and not an original, originally
47        (variables_ [i] -> Multiplicity () == 0))               // or a useless one
48      //int index = variables_ [i] -> Index ();
49      Lb (indvar) = - (Ub (indvar) = COIN_DBL_MAX);
50  }
51
52  // first initialize with values from constraints
53
54  //Jnlst()->Printf(Ipopt::J_VECTOR, J_PROBLEM, "Initial bounds for aux (initAuxs):\n");
55
56  for (std::vector <CouenneConstraint *>::const_iterator con = constraints_.begin ();
57       con != constraints_.end (); ++con) {
58
59    CouNumber
60      lb = (*((*con) -> Lb ())) (),
61      ub = (*((*con) -> Ub ())) ();
62
63    int index = (*con) -> Body () -> Index ();
64
65    assert (index >= 0);
66
67    if ((Lb (index) = CoinMax (Lb (index), lb)) <= -COUENNE_INFINITY) Lb (index) = -COIN_DBL_MAX;
68    if ((Ub (index) = CoinMin (Ub (index), ub)) >=  COUENNE_INFINITY) Ub (index) =  COIN_DBL_MAX;
69  }
70
71  // only one loop is sufficient here, since auxiliary variable are
72  // defined in such a way that w_i does NOT depend on w_j if i<j.
73
74  for (int j=0, i=nVars (); i--; j++) {
75
76    int ord = numbering_ [j];
77
78    // ignore these variables!
79    if (variables_ [ord] -> Multiplicity () == 0) {
80      Lb (ord) = - (Ub (ord) = COIN_DBL_MAX);
81      X (ord) = 0.;
82      continue;
83    }
84
85    // and handle only those with nonzero multiplicity
86    if (variables_ [ord] -> Type () == AUX) {
87
88      Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, 
89                          "w_%04d [%10g,%10g] ", ord, Lb (ord), Ub (ord));
90
91      CouNumber l, u;
92
93      variables_ [ord] -> Image () -> getBounds (l, u);
94
95      /*printf ("printing bounds: [%g %g]\n", Lb (ord), Ub (ord));
96      variables_ [ord] -> Lb () -> print (); printf ("\n");
97      variables_ [ord] -> Ub () -> print (); printf ("\n");*/
98
99      Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, 
100                          " [ --> w_%04d [%10g,%10g] ] vs [%10g %10g]", 
101                          ord, l, u, Lb (ord), Ub (ord));
102
103      // set bounds
104      if ((Lb (ord) = CoinMax (Lb (ord), l)) <= -COUENNE_INFINITY) Lb (ord) = -COIN_DBL_MAX;
105      if ((Ub (ord) = CoinMin (Ub (ord), u)) >=  COUENNE_INFINITY) Ub (ord) =  COIN_DBL_MAX;
106      //if ((lb_ [ord] = (*(aux -> Lb ())) ()) <= -COUENNE_INFINITY) lb_ [ord] = -DBL_MAX;
107      //if ((ub_ [ord] = (*(aux -> Ub ())) ()) >=  COUENNE_INFINITY) ub_ [ord] =  DBL_MAX;
108
109      Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, 
110                          " --> [%10g,%10g]\n", Lb (ord), Ub (ord));
111
112      bool integer = variables_ [ord] -> isInteger ();
113
114      if (integer) {
115        Lb (ord) = ceil  (Lb (ord) - COUENNE_EPS);
116        Ub (ord) = floor (Ub (ord) + COUENNE_EPS);
117      }
118
119      X (ord) = CoinMax (Lb (ord), CoinMin (Ub (ord), (*(variables_ [ord] -> Image ())) ()));
120    }
121  }
122}
123
124
125/// get auxiliary variables from original variables in the nonlinear
126/// problem
127void CouenneProblem::getAuxs (CouNumber * x) const {
128
129  // set point at x, don't copy
130  domain_.push (nVars (), x, domain_.lb (), domain_.ub (), false);
131
132  // set auxiliary w to f(x). This procedure is exact even though the
133  // auxiliary variables have an incomplete image, i.e. they have been
134  // decomposed previously, since they are updated with increasing
135  // index.
136  for (int j=0, i=nVars (); i--; j++) {
137
138    int index = numbering_ [j];
139    exprVar *var = variables_ [index];
140
141    if (var -> Multiplicity () > 0) {
142
143      CouNumber l, u;
144
145      if (var -> Type () == AUX)
146        var -> Image () -> getBounds (l,u);
147      else {
148        l = Lb (index);
149        u = Ub (index);
150      }
151
152      if (var -> Type () == AUX)
153        //x [index] = //not necessary, addresses of x and X are equal
154        X (index) = 
155          CoinMax (l, CoinMin (u, (*(var -> Image ())) ()));
156    } else X (index) = 0.;
157  }
158
159  domain_.pop ();
160}
161
162
163/// fill obj vector with coefficient of the (linearized) obj function
164/// (depends on sense of optimization -- invert if sense()==MAXIMIZE)
165
166void CouenneProblem::fillObjCoeff (double *&obj) {
167
168  // linearized objective can be an exprAux, an exprSub, an exprGroup,
169  // or an exprSum. In the last two cases, the components are
170  // variables or constants
171
172  expression *body = objectives_ [0] -> Body ();
173  //int sense = objectives_ [0] -> Sense ();
174
175  switch (body -> code ()) {
176
177  case COU_EXPRVAR:   //
178    obj [body -> Index ()] = 1; //(sense == MINIMIZE) ? 1 : -1;
179    break;
180
181  case COU_EXPRSUB: { //
182
183    expression **arglist = body -> ArgList ();
184
185    obj [arglist [0] -> Index ()] =  1; //(sense == MINIMIZE) ?  1 : -1;
186    obj [arglist [1] -> Index ()] = -1; //(sense == MINIMIZE) ? -1 :  1;
187
188  } break;
189
190  case COU_EXPRGROUP: { //
191
192    exprGroup *eg    = dynamic_cast <exprGroup *> (body);
193
194    const exprGroup::lincoeff &lcoe = eg -> lcoeff ();
195
196    //    if (sense == MINIMIZE) while (*index >= 0) obj [*index++] =  *coeff++;
197    //    else                   while (*index >= 0) obj [*index++] = -*coeff++;     
198
199    for (int n = lcoe.size (), i=0; n--; i++)
200      //exprGroup::lincoeff::iterator el = lcoe.begin (); el != lcoe.end (); ++el)
201      obj [lcoe [i]. first -> Index ()] = lcoe [i]. second;
202    //(sense == MINIMIZE) ?
203    //(lcoe [i]. second) :
204    //-(lcoe [i]. second);
205
206  } // no break, as exprGroup is derived from exprSum
207
208  case COU_EXPRSUM: { //
209
210    expression **arglist = body -> ArgList ();
211
212    for (int i = body -> nArgs (); i--;)
213      switch ((arglist [i]) -> code ()) {
214
215      case COU_EXPRCONST: 
216        break;
217
218      case COU_EXPRVAR: 
219        obj [arglist [i] -> Index ()] = 1; //(sense == MINIMIZE) ? 1 : -1;
220        break;
221
222      case COU_EXPRMUL: {
223
224        expression **mulArgList = arglist [i] -> ArgList ();
225        int index = mulArgList [0] -> Index ();
226
227        if (index >= 0) obj [index]                      = mulArgList [1] -> Value ();
228        else            obj [mulArgList [1] -> Index ()] = mulArgList [0] -> Value ();
229      } break;
230
231      default: 
232        Jnlst()->Printf(Ipopt::J_ERROR, J_PROBLEM,
233                        "Couenne: invalid element of sum\nAborting\n");
234        exit (-1);
235      }
236  } break;
237
238  case COU_EXPRCONST: break; // a constant objective
239
240  default:
241    Jnlst()->Printf(Ipopt::J_WARNING, J_PROBLEM,
242                    "Couenne: warning, objective function not recognized\n");
243    break;
244  }
245}
246
247
248/// set cutoff from NLP solution
249void CouenneProblem::setCutOff (CouNumber cutoff) const {
250
251  int indobj = objectives_ [0] -> Body () -> Index ();
252
253  // AW: Should we use the value of the objective variable computed by
254  //     Couenne here?
255  if ((indobj >= 0) && (cutoff < pcutoff_ -> getCutOff () - COUENNE_EPS)) {
256
257    Jnlst () -> Printf (Ipopt::J_DETAILED, J_PROBLEM,
258                        "Setting new cutoff %.10e for optimization variable index %d val = %.10e\n",
259                        cutoff, indobj,
260                        pcutoff_ -> getCutOff ());
261
262    if (Var (indobj) -> isInteger ())
263      pcutoff_    -> setCutOff (floor (cutoff + COUENNE_EPS));
264    else pcutoff_ -> setCutOff (cutoff + SafeCutoff * fabs (1. + cutoff));
265  }
266} // tolerance needed to retain feasibility
267
268
269/// Tell problem that auxiliary related to obj has a cutoff, to be
270/// used in bound tightening
271void CouenneProblem::installCutOff () const {
272
273  int indobj = objectives_ [0] -> Body () -> Index ();
274
275  if (indobj >= 0) {
276
277    // all problem are assumed to be minimization
278    double cutoff = pcutoff_ -> getCutOff();
279
280    if (cutoff < Ub (indobj))
281      Ub (indobj) = cutoff;
282
283  } else jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, 
284                           "Warning, could not install cutoff - negative objective index\n");
285}
286
287
288// clear all spurious variables pointers not referring to the variables_ vector
289void CouenneProblem::realign () {
290
291  // link variables to problem's domain
292  for (std::vector <exprVar *>::iterator i = variables_.begin ();
293       i != variables_.end (); ++i) {
294
295    (*i) -> linkDomain (&domain_);
296    (*i) -> realign (this);
297    if ((*i) -> Type () == AUX)
298      (*i) -> Image () -> realign (this);
299  }
300
301  // link variables to problem's domain
302  for (std::vector <CouenneObjective *>::iterator i = objectives_.begin ();
303       i != objectives_.end (); ++i) 
304    (*i) -> Body () -> realign (this);
305
306
307  // link variables to problem's domain
308  for (std::vector <CouenneConstraint *>::iterator i = constraints_.begin ();
309       i != constraints_.end (); ++i)
310    (*i) -> Body () -> realign (this);
311}
312
313
314/// Add list of options to be read from file
315void CouenneProblem::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
316
317  roptions -> SetRegisteringCategory ("Couenne options", Bonmin::RegisteredOptions::CouenneCategory);
318
319  roptions -> AddNumberOption
320    ("art_cutoff",
321     "Artificial cutoff",
322     COIN_DBL_MAX,
323     "Default value is infinity.");
324
325  roptions -> AddNumberOption
326    ("opt_window",
327     "Window around known optimum",
328     COIN_DBL_MAX,
329     "Default value is infinity.");
330
331  roptions -> AddNumberOption
332    ("feas_tolerance",
333     "Tolerance for constraints/auxiliary variables",
334     feas_tolerance_default,
335     "Default value is zero.");
336
337  roptions -> AddStringOption2
338    ("feasibility_bt",
339     "Feasibility-based (cheap) bound tightening",
340     "yes",
341     "no","",
342     "yes","");
343
344  roptions -> AddStringOption2
345    ("use_quadratic",
346     "Use quadratic expressions and related exprQuad class",
347     "no",
348     "no","Use an auxiliary for each bilinear term",
349     "yes","Create one only auxiliary for a quadrati expression");
350
351  roptions -> AddStringOption2
352    ("optimality_bt",
353     "Optimality-based (expensive) bound tightening",
354     "yes",
355     "no","",
356     "yes","");
357
358  roptions -> AddLowerBoundedIntegerOption
359    ("log_num_obbt_per_level",
360     "Specify the frequency (in terms of nodes) for optimality-based bound tightening.",
361     -1,0,
362     "\
363If -1, apply at every node (expensive!). \
364If 0, apply at root node only. \
365If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
366
367  roptions -> AddStringOption2
368    ("aggressive_fbbt",
369     "Aggressive feasibility-based bound tightening (to use with NLP points)",
370     "yes",
371     "no","",
372     "yes","");
373
374  roptions -> AddLowerBoundedIntegerOption
375    ("log_num_abt_per_level",
376     "Specify the frequency (in terms of nodes) for aggressive bound tightening.",
377     -1,1,
378     "\
379If -1, apply at every node (expensive!). \
380If 0, apply at root node only. \
381If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
382
383  roptions -> AddNumberOption
384    ("art_lower",
385     "Artificial lower bound",
386     -COIN_DBL_MAX,
387     "Default value is -COIN_DBL_MAX.");
388
389  roptions -> AddStringOption3
390    ("branching_object",
391     "type of branching object for variable selection",
392     "var_obj",
393     "vt_obj",   "use Violation Transfer from Tawarmalani and Sahinidis",
394     "var_obj",  "use one object for each variable",
395     "expr_obj", "use one object for each nonlinear expression");
396}
Note: See TracBrowser for help on using the repository browser.