00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "CoinHelperFunctions.hpp"
00011
00012 #include "CouenneExprHess.hpp"
00013 #include "CouenneProblem.hpp"
00014 #include "CouenneProblemElem.hpp"
00015 #include "CouenneExprAux.hpp"
00016
00017 using namespace Couenne;
00018
00019
00020
00022 ExprHess::ExprHess ():
00023
00024 nnz_ (0),
00025 iRow_ (NULL),
00026 jCol_ (NULL),
00027 numL_ (NULL),
00028 lamI_ (NULL),
00029 expr_ (NULL) {}
00030
00031
00033 ExprHess::ExprHess (const ExprHess &rhs)
00034 {operator= (rhs);}
00035
00036
00038 ExprHess &ExprHess::operator= (const ExprHess &rhs) {
00039
00040 nnz_ = rhs. nnz_;
00041
00042 iRow_ = nnz_ && rhs.iRow_ ? (int *) malloc (nnz_ * sizeof (int)) : NULL;
00043 jCol_ = nnz_ && rhs.jCol_ ? (int *) malloc (nnz_ * sizeof (int)) : NULL;
00044 numL_ = nnz_ && rhs.numL_ ? (int *) malloc (nnz_ * sizeof (int)) : NULL;
00045
00046 CoinCopyN (rhs.iRow_, nnz_, iRow_);
00047 CoinCopyN (rhs.jCol_, nnz_, jCol_);
00048 CoinCopyN (rhs.numL_, nnz_, numL_);
00049
00050 if (nnz_) {
00051
00052 lamI_ = (int **) malloc (nnz_ * sizeof (int *));
00053 expr_ = (expression ***) malloc (nnz_ * sizeof (expression **));
00054
00055 for (int i=0; i<nnz_; i++) {
00056
00057 lamI_ [i] = CoinCopyOfArray (rhs. lamI_ [i], numL_ [i]);
00058
00059 for (int j=0; j<numL_ [i]; j++)
00060 expr_ [i] [j] = rhs. expr_ [i][j] -> clone ();
00061 }
00062 }
00063
00064 return *this;
00065 }
00066
00067
00069 ExprHess *ExprHess::clone ()
00070 {return new ExprHess (*this);}
00071
00072
00074 ExprHess::~ExprHess () {
00075
00076 if (nnz_) {
00077
00078 free (iRow_);
00079 free (jCol_);
00080
00081 for (int i=0; i<nnz_; i++) {
00082
00083 for (int j=0; j<numL_ [i]; j++)
00084 delete expr_ [i] [j];
00085
00086 free (expr_ [i]);
00087 free (lamI_ [i]);
00088 }
00089
00090 free (numL_);
00091 free (lamI_);
00092 free (expr_);
00093 }
00094 }
00095
00096
00098 void HessElemFill (int i, int level,
00099 std::set <int> &list,
00100 expression *expr,
00101 int *row, int **lam, expression ***eee,
00102 CouenneProblem *,
00103 std::set <int> &globList);
00104
00105
00107 static void reAlloc (int nCur, int &nMax, int *&r, int *&c, int *&n, int **&l, expression ***&e);
00108
00109
00111 ExprHess::ExprHess (CouenneProblem *p):
00112
00113 nnz_ (0),
00114 iRow_ (NULL),
00115 jCol_ (NULL),
00116 numL_ (NULL),
00117 lamI_ (NULL),
00118 expr_ (NULL) {
00119
00120 #ifdef DEBUG
00121 printf ("creating Hessian\n");
00122
00123 #endif
00124
00127
00128 std::set <int> *deplist = new std::set <int> [1 + p -> nVars () + p -> nCons ()];
00129
00130 int nLevels = 0;
00131
00132
00133
00134
00135
00136
00137 p -> Obj (0) -> Body () -> DepList (deplist [nLevels++], STOP_AT_AUX);
00138
00139
00140
00141 for (int i = 0; i < p -> nCons (); i++) {
00142
00143 expression *c = p -> Con (i) -> Body ();
00144
00145 enum nodeType ntype = c -> Type ();
00146
00147 if (ntype == AUX ||
00148 ntype == VAR ||
00149 ntype == CONST)
00150 continue;
00151
00152 c -> DepList (deplist [nLevels++], STOP_AT_AUX);
00153 }
00154
00155
00156
00157 for (int i = 0; i < p -> nVars (); i++) {
00158
00159 exprVar *e = p -> Var (i);
00160
00161 if ((e -> Type () != AUX) ||
00162 (e -> Multiplicity () <= 0))
00163 continue;
00164
00165 e -> Image () -> DepList (deplist [nLevels++], STOP_AT_AUX);
00166 }
00167
00168 #ifdef DEBUG
00169 for (int i=0; i<nLevels; i++) {
00170 printf ("level %d: ", i);
00171 for (std::set <int>::iterator j = deplist [i]. begin () ; j != deplist [i].end (); ++j)
00172 printf ("%d ", *j);
00173 printf ("\n");
00174 }
00175 #endif
00176
00184
00185 int
00186 nVars = p -> nVars (),
00187 curSize = 0;
00188
00190 for (int i=0; i < nVars; i++) {
00191
00192
00193 int *rnz = (int *) malloc (nVars * sizeof (int));
00194 int **lam = (int **) malloc (nVars * sizeof (int *));
00195 expression ***eee = (expression ***) malloc (nVars * sizeof (expression **));
00196
00197 std::set <int> globDepList;
00198
00199 CoinFillN (rnz, nVars, 0);
00200
00201
00202 for (int j=nVars; j--;) {
00203 *lam++ = NULL;
00204 *eee++ = NULL;
00205 }
00206
00207 lam -= nVars;
00208 eee -= nVars;
00209
00210
00211 int level = 0;
00212
00214 HessElemFill (i, 0, deplist [0], p -> Obj (0) -> Body (), rnz, lam, eee, p, globDepList);
00215
00216 ++level;
00217
00218 for (int j = 0; j < p -> nCons (); j++) {
00219
00220 CouenneConstraint *c = p -> Con (j);
00221 enum nodeType ctype = c -> Body () -> Type ();
00222
00223 if (ctype == AUX ||
00224 ctype == VAR)
00225 continue;
00226
00227 HessElemFill (i, level, deplist [level], c -> Body (), rnz, lam, eee, p, globDepList);
00228
00229 ++level;
00230 }
00231
00232
00233
00234 for (int j = 0; j < p -> nVars (); j++) {
00235
00236 exprVar *e = p -> Var (j);
00237
00238 if ((e -> Type () != AUX) ||
00239 (e -> Multiplicity () <= 0))
00240 continue;
00241
00242 HessElemFill (i, level, deplist [level], e -> Image (), rnz, lam, eee, p, globDepList);
00243
00244 ++level;
00245 }
00246
00247
00248
00249 for (std::set <int>::iterator j = globDepList.begin (); j != globDepList. end (); ++j) {
00250
00251 assert (*j <= i);
00252 assert (rnz [*j]);
00253
00254 reAlloc (nnz_ + 1, curSize, iRow_, jCol_, numL_, lamI_, expr_);
00255
00256 iRow_ [nnz_] = i;
00257 jCol_ [nnz_] = *j;
00258 numL_ [nnz_] = rnz [*j];
00259 lamI_ [nnz_] = (int *) realloc (lam [*j], rnz [*j] * sizeof (int));
00260 expr_ [nnz_] = (expression **) realloc (eee [*j], rnz [*j] * sizeof (expression *));
00261
00262 ++nnz_;
00263 }
00264
00265 free (rnz);
00266 free (lam);
00267 free (eee);
00268 }
00269
00270 delete [] deplist;
00271
00272 #ifdef DEBUG
00273 printf ("hessian: %d nonzeros\n", nnz_);
00274
00275 for (int i=0; i<nnz_; i++) {
00276
00277 printf ("[%d,%d]: %d lambdas: ",
00278 iRow_ [i], jCol_ [i], numL_ [i]);
00279
00280 fflush (stdout);
00281
00282 for (int j=0; j<numL_ [i]; j++) {
00283 printf ("(%d,", lamI_ [i][j]);
00284 fflush (stdout);
00285 expr_ [i][j] -> print ();
00286 fflush (stdout);
00287 printf (") ");
00288 }
00289 printf ("\n");
00290 }
00291 #endif
00292 }
00293
00294
00295 #define reallocStep 100
00296
00297
00298
00299 void HessElemFill (int i,
00300 int level,
00301 std::set <int> &list,
00302 expression *expr,
00303 int *nnz,
00304 int **lam,
00305 expression ***eee,
00306 CouenneProblem *p,
00307 std::set <int> &globList) {
00308
00309 if (list. find (i) == list.end () ||
00310 (expr -> Linearity () <= LINEAR))
00311 return;
00312
00313
00314
00315 expression
00316 *d1 = expr -> differentiate (i),
00317 *sd1 = d1 -> simplify (),
00318 *rd1 = (sd1 ? sd1 : d1);
00319
00320 rd1 -> realign (p);
00321
00322 for (std::set <int>::iterator k = list.begin (); k != list. end (); ++k) {
00323
00324 if (*k > i)
00325 continue;
00326
00327
00328
00329 expression
00330 *d2 = rd1 -> differentiate (*k),
00331 *sd2 = d2 -> simplify (),
00332 *rd2 = (sd2 ? sd2 : d2);
00333
00334 #ifdef DEBUG
00335 printf (" rd2 [x_%d, x_%d]: d/d x_%d = ", *k, i, *k); fflush (stdout);
00336 rd1 -> print (); printf (" -> d/(d x_%d,d x_%d) = ", *k, i);
00337 rd2 -> print (); printf ("\n");
00338 #endif
00339
00340 if (sd2) delete d2;
00341
00342 if ((rd2 -> Type () != CONST) ||
00343 (rd2 -> Value () != 0.)) {
00344
00345
00346
00347 int &curNNZ = nnz [*k];
00348
00349 if (!curNNZ &&
00350 globList.find (*k) == globList. end ())
00351 globList.insert (*k);
00352
00353 if (!(curNNZ % reallocStep)) {
00354
00355 lam [*k] = (int *) realloc (lam [*k], (curNNZ + reallocStep) * sizeof (int));
00356 eee [*k] = (expression **) realloc (eee [*k], (curNNZ + reallocStep) * sizeof (expression *));
00357 }
00358
00359 rd2 -> realign (p);
00360
00361 lam [*k] [curNNZ] = level;
00362 eee [*k] [curNNZ++] = rd2;
00363
00364 } else delete rd2;
00365 }
00366
00367 if (sd1) delete sd1;
00368 delete d1;
00369 }
00370
00371 static void reAlloc (int nCur, int &nMax, int *&r, int *&c, int *&n, int **&l, expression ***&e) {
00372
00373 if (nCur >= nMax) {
00374
00375 nMax += reallocStep;
00376
00377 r = (int *) realloc (r, nMax * sizeof (int));
00378 c = (int *) realloc (c, nMax * sizeof (int));
00379 n = (int *) realloc (n, nMax * sizeof (int));
00380 l = (int **) realloc (l, nMax * sizeof (int *));
00381 e = (expression ***) realloc (e, nMax * sizeof (expression **));
00382 }
00383 }