1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file nlpioracle.c 26 * @ingroup OTHER_CFILES 27 * @brief implementation of NLPI oracle 28 * @author Stefan Vigerske 29 * 30 * @todo jacobi evaluation should be sparse 31 */ 32 33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 34 35 #include "scip/scip.h" 36 #include "scip/nlpioracle.h" 37 #include "scip/exprinterpret.h" 38 #include "scip/expr_pow.h" 39 #include "scip/expr_varidx.h" 40 41 #include <string.h> /* for strlen */ 42 43 /**@name NLPI Oracle data structures */ 44 /**@{ */ 45 46 struct SCIP_NlpiOracleCons 47 { 48 SCIP_Real lhs; /**< left hand side (for constraint) or constant (for objective) */ 49 SCIP_Real rhs; /**< right hand side (for constraint) or constant (for objective) */ 50 51 int linsize; /**< length of linidxs and lincoefs arrays */ 52 int nlinidxs; /**< number of linear variable indices and coefficients */ 53 int* linidxs; /**< variable indices in linear part, or NULL if none */ 54 SCIP_Real* lincoefs; /**< variable coefficients in linear part, of NULL if none */ 55 56 SCIP_EXPR* expr; /**< expression for nonlinear part, or NULL if none */ 57 SCIP_EXPRINTDATA* exprintdata; /**< expression interpret data for expression, or NULL if no expr or not compiled yet */ 58 59 char* name; /**< name of constraint */ 60 }; 61 typedef struct SCIP_NlpiOracleCons SCIP_NLPIORACLECONS; 62 63 struct SCIP_NlpiOracle 64 { 65 char* name; /**< name of problem */ 66 67 int varssize; /**< length of variables related arrays */ 68 int nvars; /**< number of variables */ 69 SCIP_Real* varlbs; /**< array with variable lower bounds */ 70 SCIP_Real* varubs; /**< array with variable upper bounds */ 71 char** varnames; /**< array with variable names */ 72 int* varlincount; /**< array with number of appearances of variable in linear part of objective or constraints */ 73 int* varnlcount; /**< array with number of appearances of variable in nonlinear part of objective or constraints */ 74 75 int consssize; /**< length of constraints related arrays */ 76 int nconss; /**< number of constraints */ 77 SCIP_NLPIORACLECONS** conss; /**< constraints, or NULL if none */ 78 79 SCIP_NLPIORACLECONS* objective; /**< objective */ 80 81 int* jacoffsets; /**< rowwise jacobi sparsity pattern: constraint offsets in jaccols */ 82 int* jaccols; /**< rowwise jacobi sparsity pattern: indices of variables appearing in constraints */ 83 84 int* heslagoffsets; /**< rowwise sparsity pattern of hessian matrix of Lagrangian: row offsets in heslagcol */ 85 int* heslagcols; /**< rowwise sparsity pattern of hessian matrix of Lagrangian: column indices; sorted for each row */ 86 87 SCIP_EXPRINT* exprinterpreter; /**< interpreter for expressions: evaluation and derivatives */ 88 SCIP_CLOCK* evalclock; /**< clock measuring evaluation time */ 89 }; 90 91 /**@} */ 92 93 /*lint -e440*/ 94 /*lint -e441*/ 95 /*lint -e866*/ 96 97 /**@name Local functions */ 98 /**@{ */ 99 100 /** ensures that those arrays in oracle that store information on variables have at least a given length */ 101 static 102 SCIP_RETCODE ensureVarsSize( 103 SCIP* scip, /**< SCIP data structure */ 104 SCIP_NLPIORACLE* oracle, /**< NLPIORACLE data structure */ 105 int minsize /**< minimal required size */ 106 ) 107 { 108 assert(oracle != NULL); 109 110 if( minsize > oracle->varssize ) 111 { 112 int newsize; 113 114 newsize = SCIPcalcMemGrowSize(scip, minsize); 115 assert(newsize >= minsize); 116 117 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &oracle->varlbs, oracle->varssize, newsize) ); 118 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &oracle->varubs, oracle->varssize, newsize) ); 119 if( oracle->varnames != NULL ) 120 { 121 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &oracle->varnames, oracle->varssize, newsize) ); 122 } 123 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &oracle->varlincount, oracle->varssize, newsize) ); 124 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &oracle->varnlcount, oracle->varssize, newsize) ); 125 126 oracle->varssize = newsize; 127 } 128 assert(oracle->varssize >= minsize); 129 130 return SCIP_OKAY; 131 } 132 133 /** ensures that constraints array in oracle has at least a given length */ 134 static 135 SCIP_RETCODE ensureConssSize( 136 SCIP* scip, /**< SCIP data structure */ 137 SCIP_NLPIORACLE* oracle, /**< NLPIORACLE data structure */ 138 int minsize /**< minimal required size */ 139 ) 140 { 141 assert(oracle != NULL); 142 143 SCIP_CALL( SCIPensureBlockMemoryArray(scip, &oracle->conss, &oracle->consssize, minsize) ); 144 assert(oracle->consssize >= minsize); 145 146 return SCIP_OKAY; 147 } 148 149 /** ensures that arrays for linear part in a oracle constraints have at least a given length */ 150 static 151 SCIP_RETCODE ensureConsLinSize( 152 SCIP* scip, /**< SCIP data structure */ 153 SCIP_NLPIORACLECONS* cons, /**< oracle constraint */ 154 int minsize /**< minimal required size */ 155 ) 156 { 157 assert(cons != NULL); 158 159 if( minsize > cons->linsize ) 160 { 161 int newsize; 162 163 newsize = SCIPcalcMemGrowSize(scip, minsize); 164 assert(newsize >= minsize); 165 166 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &cons->linidxs, cons->linsize, newsize) ); 167 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &cons->lincoefs, cons->linsize, newsize) ); 168 cons->linsize = newsize; 169 } 170 assert(cons->linsize >= minsize); 171 172 return SCIP_OKAY; 173 } 174 175 /** ensures that a given array of integers has at least a given length */ 176 static 177 SCIP_RETCODE ensureIntArraySize( 178 SCIP* scip, /**< SCIP data structure */ 179 int** intarray, /**< array of integers */ 180 int* len, /**< length of array (modified if reallocated) */ 181 int minsize /**< minimal required array length */ 182 ) 183 { 184 assert(intarray != NULL); 185 assert(len != NULL); 186 187 SCIP_CALL( SCIPensureBlockMemoryArray(scip, intarray, len, minsize) ); 188 assert(*len >= minsize); 189 190 return SCIP_OKAY; 191 } 192 193 /** Invalidates the sparsity pattern of the Jacobian. 194 * Should be called when constraints are added or deleted. 195 */ 196 static 197 void invalidateJacobiSparsity( 198 SCIP* scip, /**< SCIP data structure */ 199 SCIP_NLPIORACLE* oracle /**< pointer to store NLPIORACLE data structure */ 200 ) 201 { 202 assert(oracle != NULL); 203 204 SCIPdebugMessage("%p invalidate jacobian sparsity\n", (void*)oracle); 205 206 if( oracle->jacoffsets == NULL ) 207 { /* nothing to do */ 208 assert(oracle->jaccols == NULL); 209 return; 210 } 211 212 assert(oracle->jaccols != NULL); 213 SCIPfreeBlockMemoryArray(scip, &oracle->jaccols, oracle->jacoffsets[oracle->nconss]); 214 SCIPfreeBlockMemoryArray(scip, &oracle->jacoffsets, oracle->nconss + 1); 215 } 216 217 /** Invalidates the sparsity pattern of the Hessian of the Lagragian. 218 * Should be called when the objective is set or constraints are added or deleted. 219 */ 220 static 221 void invalidateHessianLagSparsity( 222 SCIP* scip, /**< SCIP data structure */ 223 SCIP_NLPIORACLE* oracle /**< pointer to store NLPIORACLE data structure */ 224 ) 225 { 226 assert(oracle != NULL); 227 228 SCIPdebugMessage("%p invalidate hessian lag sparsity\n", (void*)oracle); 229 230 if( oracle->heslagoffsets == NULL ) 231 { /* nothing to do */ 232 assert(oracle->heslagcols == NULL); 233 return; 234 } 235 236 assert(oracle->heslagcols != NULL); 237 SCIPfreeBlockMemoryArray(scip, &oracle->heslagcols, oracle->heslagoffsets[oracle->nvars]); 238 SCIPfreeBlockMemoryArray(scip, &oracle->heslagoffsets, oracle->nvars + 1); 239 } 240 241 /** increases or decreases variable counts in oracle w.r.t. linear and nonlinear appearance */ 242 static 243 SCIP_RETCODE updateVariableCounts( 244 SCIP* scip, /**< SCIP data structure */ 245 SCIP_NLPIORACLE* oracle, /**< oracle data structure */ 246 int factor, /**< whether to add (factor=1) or remove (factor=-1) variable counts */ 247 int nlinidxs, /**< number of linear indices */ 248 const int* linidxs, /**< indices of variables in linear part */ 249 SCIP_EXPR* expr /**< expression */ 250 ) 251 { 252 int j; 253 254 assert(oracle != NULL); 255 assert(oracle->varlincount != NULL || (nlinidxs == 0 && expr == NULL)); 256 assert(oracle->varnlcount != NULL || (nlinidxs == 0 && expr == NULL)); 257 assert(factor == 1 || factor == -1); 258 assert(nlinidxs == 0 || linidxs != NULL); 259 260 for( j = 0; j < nlinidxs; ++j ) 261 { 262 oracle->varlincount[linidxs[j]] += factor; 263 assert(oracle->varlincount[linidxs[j]] >= 0); 264 } 265 266 if( expr != NULL ) 267 { 268 SCIP_EXPRITER* it; 269 270 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 271 SCIP_CALL( SCIPexpriterInit(it, expr, SCIP_EXPRITER_DFS, FALSE) ); 272 273 for( ; !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) 274 if( SCIPisExprVaridx(scip, expr) ) 275 { 276 oracle->varnlcount[SCIPgetIndexExprVaridx(expr)] += factor; 277 assert(oracle->varnlcount[SCIPgetIndexExprVaridx(expr)] >= 0); 278 } 279 280 SCIPfreeExpriter(&it); 281 } 282 283 return SCIP_OKAY; 284 } 285 286 /** sorts a linear term, merges duplicate entries and removes entries with coefficient 0.0 */ 287 static 288 void sortLinearCoefficients( 289 int* nidxs, /**< number of variables */ 290 int* idxs, /**< indices of variables */ 291 SCIP_Real* coefs /**< coefficients of variables */ 292 ) 293 { 294 int offset; 295 int j; 296 297 assert(nidxs != NULL); 298 assert(idxs != NULL || *nidxs == 0); 299 assert(coefs != NULL || *nidxs == 0); 300 301 if( *nidxs == 0 ) 302 return; 303 304 SCIPsortIntReal(idxs, coefs, *nidxs); 305 306 offset = 0; 307 j = 0; 308 while( j+offset < *nidxs ) 309 { 310 assert(idxs[j] >= 0); /*lint !e613*/ 311 312 /* move j+offset to j, if different */ 313 if( offset > 0 ) 314 { 315 idxs[j] = idxs[j+offset]; /*lint !e613*/ 316 coefs[j] = coefs[j+offset]; /*lint !e613*/ 317 } 318 319 /* add up coefs for j+offset+1... as long as they have the same index */ 320 while( j+offset+1 < *nidxs && idxs[j] == idxs[j+offset+1] ) /*lint !e613*/ 321 { 322 coefs[j] += coefs[j+offset+1]; /*lint !e613*/ 323 ++offset; 324 } 325 326 /* if j'th element is 0, increase offset, otherwise increase j */ 327 if( coefs[j] == 0.0 ) /*lint !e613*/ 328 ++offset; 329 else 330 ++j; 331 } 332 *nidxs -= offset; 333 } 334 335 /** creates a NLPI constraint from given constraint data */ 336 static 337 SCIP_RETCODE createConstraint( 338 SCIP* scip, /**< SCIP data structure */ 339 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 340 SCIP_NLPIORACLECONS** cons, /**< buffer where to store pointer to constraint */ 341 int nlinidxs, /**< length of linear part */ 342 const int* linidxs, /**< indices of linear part, or NULL if nlinidxs == 0 */ 343 const SCIP_Real* lincoefs, /**< coefficients of linear part, or NULL if nlinidxs == 0 */ 344 SCIP_EXPR* expr, /**< expression, or NULL */ 345 SCIP_Real lhs, /**< left-hand-side of constraint */ 346 SCIP_Real rhs, /**< right-hand-side of constraint */ 347 const char* name /**< name of constraint, or NULL */ 348 ) 349 { 350 assert(cons != NULL); 351 assert(nlinidxs >= 0); 352 assert(linidxs != NULL || nlinidxs == 0); 353 assert(lincoefs != NULL || nlinidxs == 0); 354 assert(EPSLE(lhs, rhs, SCIP_DEFAULT_EPSILON)); 355 356 SCIP_CALL( SCIPallocClearBlockMemory(scip, cons) ); 357 assert(*cons != NULL); 358 359 if( nlinidxs > 0 ) 360 { 361 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*cons)->linidxs, linidxs, nlinidxs) ); 362 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*cons)->lincoefs, lincoefs, nlinidxs) ); 363 (*cons)->linsize = nlinidxs; 364 (*cons)->nlinidxs = nlinidxs; 365 366 /* sort, merge duplicates, remove zero's */ 367 sortLinearCoefficients(&(*cons)->nlinidxs, (*cons)->linidxs, (*cons)->lincoefs); 368 assert((*cons)->linidxs[0] >= 0); 369 } 370 371 if( expr != NULL ) 372 { 373 (*cons)->expr = expr; 374 SCIPcaptureExpr(expr); 375 376 SCIP_CALL( SCIPexprintCompile(scip, oracle->exprinterpreter, (*cons)->expr, &(*cons)->exprintdata) ); 377 } 378 379 if( lhs > rhs ) 380 { 381 assert(EPSEQ(lhs, rhs, SCIP_DEFAULT_EPSILON)); 382 lhs = rhs; 383 } 384 (*cons)->lhs = lhs; 385 (*cons)->rhs = rhs; 386 387 if( name != NULL ) 388 { 389 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*cons)->name, name, strlen(name)+1) ); 390 } 391 392 /* add variable counts */ 393 SCIP_CALL( updateVariableCounts(scip, oracle, 1, (*cons)->nlinidxs, (*cons)->linidxs, (*cons)->expr) ); 394 395 return SCIP_OKAY; 396 } 397 398 /** frees a constraint */ 399 static 400 SCIP_RETCODE freeConstraint( 401 SCIP* scip, /**< SCIP data structure */ 402 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 403 SCIP_NLPIORACLECONS** cons, /**< pointer to constraint that should be freed */ 404 SCIP_Bool updatevarcount /**< whether the update variable counts (typically TRUE) */ 405 ) 406 { 407 assert(oracle != NULL); 408 assert(cons != NULL); 409 assert(*cons != NULL); 410 411 SCIPdebugMessage("free constraint %p\n", (void*)*cons); 412 413 /* remove variable counts */ 414 if( updatevarcount ) 415 { 416 SCIP_CALL( updateVariableCounts(scip, oracle, -1, (*cons)->nlinidxs, (*cons)->linidxs, (*cons)->expr) ); 417 } 418 419 SCIPfreeBlockMemoryArrayNull(scip, &(*cons)->linidxs, (*cons)->linsize); 420 SCIPfreeBlockMemoryArrayNull(scip, &(*cons)->lincoefs, (*cons)->linsize); 421 422 if( (*cons)->expr != NULL ) 423 { 424 SCIP_CALL( SCIPexprintFreeData(scip, oracle->exprinterpreter, (*cons)->expr, &(*cons)->exprintdata) ); 425 SCIP_CALL( SCIPreleaseExpr(scip, &(*cons)->expr) ); 426 } 427 428 if( (*cons)->name != NULL ) 429 { 430 SCIPfreeBlockMemoryArrayNull(scip, &(*cons)->name, strlen((*cons)->name)+1); 431 } 432 433 SCIPfreeBlockMemory(scip, cons); 434 assert(*cons == NULL); 435 436 return SCIP_OKAY; 437 } 438 439 /** frees all constraints 440 * 441 * \attention This omits updating the variable counts in the oracle. 442 */ 443 static 444 SCIP_RETCODE freeConstraints( 445 SCIP* scip, /**< SCIP data structure */ 446 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 447 ) 448 { 449 int i; 450 451 assert(oracle != NULL); 452 453 SCIPdebugMessage("%p free constraints\n", (void*)oracle); 454 455 for( i = 0; i < oracle->nconss; ++i ) 456 { 457 SCIP_CALL( freeConstraint(scip, oracle, &oracle->conss[i], FALSE) ); 458 assert(oracle->conss[i] == NULL); 459 } 460 oracle->nconss = 0; 461 462 SCIPfreeBlockMemoryArrayNull(scip, &oracle->conss, oracle->consssize); 463 oracle->consssize = 0; 464 465 return SCIP_OKAY; 466 } 467 468 /** moves one variable 469 * The place where it moves to need to be empty (all NULL) but allocated. 470 * Note that this function does not update the variable indices in the constraints! 471 */ 472 static 473 SCIP_RETCODE moveVariable( 474 SCIP* scip, /**< SCIP data structure */ 475 SCIP_NLPIORACLE* oracle, /**< pointer to store NLPIORACLE data structure */ 476 int fromidx, /**< index of variable to move */ 477 int toidx /**< index of place where to move variable to */ 478 ) 479 { 480 assert(oracle != NULL); 481 482 SCIPdebugMessage("%p move variable\n", (void*)oracle); 483 484 assert(0 <= fromidx); 485 assert(0 <= toidx); 486 assert(fromidx < oracle->nvars); 487 assert(toidx < oracle->nvars); 488 489 assert(oracle->varnames == NULL || oracle->varnames[toidx] == NULL); 490 491 oracle->varlbs[toidx] = oracle->varlbs[fromidx]; 492 oracle->varubs[toidx] = oracle->varubs[fromidx]; 493 oracle->varlbs[fromidx] = -SCIPinfinity(scip); 494 oracle->varubs[fromidx] = SCIPinfinity(scip); 495 496 oracle->varlincount[toidx] = oracle->varlincount[fromidx]; 497 oracle->varnlcount[toidx] = oracle->varnlcount[fromidx]; 498 oracle->varlincount[fromidx] = 0; 499 oracle->varnlcount[fromidx] = 0; 500 501 if( oracle->varnames != NULL ) 502 { 503 oracle->varnames[toidx] = oracle->varnames[fromidx]; 504 oracle->varnames[fromidx] = NULL; 505 } 506 507 return SCIP_OKAY; 508 } 509 510 /** frees all variables */ 511 static 512 void freeVariables( 513 SCIP* scip, /**< SCIP data structure */ 514 SCIP_NLPIORACLE* oracle /**< pointer to store NLPIORACLE data structure */ 515 ) 516 { 517 int i; 518 519 assert(oracle != NULL); 520 521 SCIPdebugMessage("%p free variables\n", (void*)oracle); 522 523 if( oracle->varnames != NULL ) 524 { 525 for( i = 0; i < oracle->nvars; ++i ) 526 { 527 if( oracle->varnames[i] != NULL ) 528 { 529 SCIPfreeBlockMemoryArray(scip, &oracle->varnames[i], strlen(oracle->varnames[i])+1); /*lint !e866*/ 530 } 531 } 532 SCIPfreeBlockMemoryArrayNull(scip, &oracle->varnames, oracle->varssize); 533 } 534 oracle->nvars = 0; 535 536 SCIPfreeBlockMemoryArrayNull(scip, &oracle->varlbs, oracle->varssize); 537 SCIPfreeBlockMemoryArrayNull(scip, &oracle->varubs, oracle->varssize); 538 SCIPfreeBlockMemoryArrayNull(scip, &oracle->varlincount, oracle->varssize); 539 SCIPfreeBlockMemoryArrayNull(scip, &oracle->varnlcount, oracle->varssize); 540 541 oracle->varssize = 0; 542 } 543 544 /** applies a mapping of indices to one array of indices */ 545 static 546 void mapIndices( 547 int* indexmap, /**< mapping from old variable indices to new indices */ 548 int nindices, /**< number of indices in indices1 and indices2 */ 549 int* indices /**< array of indices to adjust */ 550 ) 551 { 552 assert(indexmap != NULL); 553 assert(nindices == 0 || indices != NULL); 554 555 for( ; nindices ; --nindices, ++indices ) 556 { 557 assert(indexmap[*indices] >= 0); 558 *indices = indexmap[*indices]; 559 } 560 } 561 562 /** removes entries with index -1 (marked as deleted) from array of linear elements 563 * assumes that array is sorted by index, i.e., all -1 are at the beginning 564 */ 565 static 566 void clearDeletedLinearElements( 567 int** linidxs, /**< variable indices */ 568 SCIP_Real** coefs, /**< variable coefficients */ 569 int* nidxs /**< number of indices */ 570 ) 571 { 572 int i; 573 int offset; 574 575 SCIPdebugMessage("clear deleted linear elements\n"); 576 577 assert(linidxs != NULL); 578 assert(*linidxs != NULL); 579 assert(coefs != NULL); 580 assert(*coefs != NULL); 581 assert(nidxs != NULL); 582 assert(*nidxs > 0); 583 584 /* search for beginning of non-delete entries @todo binary search? */ 585 for( offset = 0; offset < *nidxs; ++offset ) 586 if( (*linidxs)[offset] >= 0 ) 587 break; 588 589 /* nothing was deleted */ 590 if( offset == 0 ) 591 return; 592 593 /* some or all elements were deleted -> move remaining ones front */ 594 for( i = 0; i < *nidxs - offset; ++i ) 595 { 596 (*linidxs)[i] = (*linidxs)[i+offset]; 597 (*coefs)[i] = (*coefs) [i+offset]; 598 } 599 *nidxs -= offset; 600 } 601 602 /** computes the value of a function */ 603 static 604 SCIP_RETCODE evalFunctionValue( 605 SCIP* scip, /**< SCIP data structure */ 606 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 607 SCIP_NLPIORACLECONS* cons, /**< oracle constraint */ 608 const SCIP_Real* x, /**< the point where to evaluate */ 609 SCIP_Real* val /**< pointer to store function value */ 610 ) 611 { /*lint --e{715}*/ 612 assert(oracle != NULL); 613 assert(cons != NULL); 614 assert(x != NULL || oracle->nvars == 0); 615 assert(val != NULL); 616 617 SCIPdebugMessage("%p eval function value\n", (void*)oracle); 618 619 *val = 0.0; 620 621 if( cons->nlinidxs > 0 ) 622 { 623 int* linidxs; 624 SCIP_Real* lincoefs; 625 int nlin; 626 627 nlin = cons->nlinidxs; 628 linidxs = cons->linidxs; 629 lincoefs = cons->lincoefs; 630 assert(linidxs != NULL); 631 assert(lincoefs != NULL); 632 assert(x != NULL); 633 634 for( ; nlin > 0; --nlin, ++linidxs, ++lincoefs ) 635 *val += *lincoefs * x[*linidxs]; 636 } 637 638 if( cons->expr != NULL ) 639 { 640 SCIP_Real nlval; 641 642 SCIP_CALL( SCIPexprintEval(scip, oracle->exprinterpreter, cons->expr, cons->exprintdata, (SCIP_Real*)x, &nlval) ); 643 if( !SCIPisFinite(nlval) || SCIPisInfinity(scip, ABS(nlval)) ) 644 *val = nlval; 645 else 646 *val += nlval; 647 } 648 649 return SCIP_OKAY; 650 } 651 652 /** computes the value and gradient of a function 653 * 654 * @return SCIP_INVALIDDATA, if the function or its gradient could not be evaluated (domain error, etc.) 655 */ 656 static 657 SCIP_RETCODE evalFunctionGradient( 658 SCIP* scip, /**< SCIP data structure */ 659 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 660 SCIP_NLPIORACLECONS* cons, /**< oracle constraint */ 661 const SCIP_Real* x, /**< the point where to evaluate */ 662 SCIP_Bool isnewx, /**< has the point x changed since the last call to some evaluation function? */ 663 SCIP_Real* RESTRICT val, /**< pointer to store function value */ 664 SCIP_Real* RESTRICT grad /**< pointer to store function gradient */ 665 ) 666 { /*lint --e{715}*/ 667 assert(oracle != NULL); 668 assert(x != NULL || oracle->nvars == 0); 669 assert(val != NULL); 670 assert(grad != NULL); 671 672 SCIPdebugMessage("%p eval function gradient\n", (void*)oracle); 673 674 *val = 0.0; 675 BMSclearMemoryArray(grad, oracle->nvars); 676 677 if( cons->expr != NULL ) 678 { 679 SCIP_Real nlval; 680 int i; 681 682 SCIPdebugMsg(scip, "eval gradient of "); 683 SCIPdebug( if( isnewx ) {printf("\nx ="); for( i = 0; i < oracle->nvars; ++i) printf(" %g", x[i]); printf("\n");} ) 684 685 SCIP_CALL( SCIPexprintGrad(scip, oracle->exprinterpreter, cons->expr, cons->exprintdata, (SCIP_Real*)x, isnewx, &nlval, grad) ); 686 687 SCIPdebug( printf("g ="); for( i = 0; i < oracle->nvars; ++i) printf(" %g", grad[i]); printf("\n"); ) 688 689 /* check for eval error */ 690 if( !SCIPisFinite(nlval) || SCIPisInfinity(scip, ABS(nlval)) ) 691 { 692 SCIPdebugMessage("gradient evaluation yield invalid function value %g\n", nlval); 693 return SCIP_INVALIDDATA; /* indicate that the function could not be evaluated at given point */ 694 } 695 for( i = 0; i < oracle->nvars; ++i ) 696 if( !SCIPisFinite(grad[i]) ) 697 { 698 SCIPdebugMessage("gradient evaluation yield invalid gradient value %g\n", grad[i]); 699 return SCIP_INVALIDDATA; /* indicate that the function could not be evaluated at given point */ 700 } 701 702 *val += nlval; 703 } 704 705 if( cons->nlinidxs > 0 ) 706 { 707 int* linidxs; 708 SCIP_Real* lincoefs; 709 int nlin; 710 711 nlin = cons->nlinidxs; 712 linidxs = cons->linidxs; 713 lincoefs = cons->lincoefs; 714 assert(linidxs != NULL); 715 assert(lincoefs != NULL); 716 assert(x != NULL); 717 718 for( ; nlin > 0; --nlin, ++linidxs, ++lincoefs ) 719 { 720 *val += *lincoefs * x[*linidxs]; 721 grad[*linidxs] += *lincoefs; 722 } 723 } 724 725 return SCIP_OKAY; 726 } 727 728 /** collects indices of nonzero entries in the lower-left part of the hessian matrix of an expression 729 * adds the indices to a given set of indices, avoiding duplicates */ 730 static 731 SCIP_RETCODE hessLagSparsitySetNzFlagForExpr( 732 SCIP* scip, /**< SCIP data structure */ 733 SCIP_NLPIORACLE* oracle, /**< NLPI oracle */ 734 int** colnz, /**< indices of nonzero entries for each column */ 735 int* collen, /**< space allocated to store indices of nonzeros for each column */ 736 int* colnnz, /**< number of nonzero entries for each column */ 737 int* nzcount, /**< counter for total number of nonzeros; should be increased when nzflag is set to 1 the first time */ 738 SCIP_EXPR* expr, /**< expression */ 739 SCIP_EXPRINTDATA* exprintdata, /**< expression interpreter data for expression */ 740 int dim /**< dimension of matrix */ 741 ) 742 { 743 SCIP_Real* x; 744 int* rowidxs; 745 int* colidxs; 746 int nnz; 747 int row; 748 int col; 749 int pos; 750 int i; 751 752 assert(oracle != NULL); 753 assert(colnz != NULL); 754 assert(collen != NULL); 755 assert(colnnz != NULL); 756 assert(nzcount != NULL); 757 assert(expr != NULL); 758 assert(dim >= 0); 759 760 SCIPdebugMessage("%p hess lag sparsity set nzflag for expr\n", (void*)oracle); 761 762 SCIP_CALL( SCIPallocBufferArray(scip, &x, oracle->nvars) ); 763 for( i = 0; i < oracle->nvars; ++i ) 764 x[i] = 2.0; /* hope that this value does not make much trouble for the evaluation routines */ 765 766 SCIP_CALL( SCIPexprintHessianSparsity(scip, oracle->exprinterpreter, expr, exprintdata, x, &rowidxs, &colidxs, &nnz) ); 767 768 for( i = 0; i < nnz; ++i ) 769 { 770 row = rowidxs[i]; 771 col = colidxs[i]; 772 773 assert(row < oracle->nvars); 774 assert(col <= row); 775 776 if( colnz[row] == NULL || !SCIPsortedvecFindInt(colnz[row], col, colnnz[row], &pos) ) 777 { 778 SCIP_CALL( ensureIntArraySize(scip, &colnz[row], &collen[row], colnnz[row]+1) ); 779 SCIPsortedvecInsertInt(colnz[row], col, &colnnz[row], NULL); 780 ++*nzcount; 781 } 782 } 783 784 SCIPfreeBufferArray(scip, &x); 785 786 return SCIP_OKAY; 787 } 788 789 /** adds hessian of an expression into hessian structure */ 790 static 791 SCIP_RETCODE hessLagAddExpr( 792 SCIP* scip, /**< SCIP data structure */ 793 SCIP_NLPIORACLE* oracle, /**< oracle */ 794 SCIP_Real weight, /**< weight of quadratic part */ 795 const SCIP_Real* x, /**< point for which hessian should be returned */ 796 SCIP_Bool new_x, /**< whether point has been evaluated before */ 797 SCIP_EXPR* expr, /**< expression */ 798 SCIP_EXPRINTDATA* exprintdata, /**< expression interpreter data for expression */ 799 int* hesoffset, /**< row offsets in sparse matrix that is to be filled */ 800 int* hescol, /**< column indices in sparse matrix that is to be filled */ 801 SCIP_Real* values /**< buffer for values of sparse matrix that is to be filled */ 802 ) 803 { 804 SCIP_Real val; 805 SCIP_Real* h; 806 int* rowidxs; 807 int* colidxs; 808 int nnz; 809 int row; 810 int col; 811 int pos; 812 int i; 813 814 SCIPdebugMessage("%p hess lag add expr\n", (void*)oracle); 815 816 assert(oracle != NULL); 817 assert(x != NULL || new_x == FALSE); 818 assert(expr != NULL); 819 assert(hesoffset != NULL); 820 assert(hescol != NULL); 821 assert(values != NULL); 822 823 SCIP_CALL( SCIPexprintHessian(scip, oracle->exprinterpreter, expr, exprintdata, (SCIP_Real*)x, new_x, &val, &rowidxs, &colidxs, &h, &nnz) ); 824 if( !SCIPisFinite(val) ) 825 { 826 SCIPdebugMessage("hessian evaluation yield invalid function value %g\n", val); 827 return SCIP_INVALIDDATA; /* indicate that the function could not be evaluated at given point */ 828 } 829 830 for( i = 0; i < nnz; ++i ) 831 { 832 if( !SCIPisFinite(h[i]) ) 833 { 834 SCIPdebugMessage("hessian evaluation yield invalid hessian value %g\n", *h); 835 return SCIP_INVALIDDATA; /* indicate that the function could not be evaluated at given point */ 836 } 837 838 if( h[i] == 0.0 ) 839 continue; 840 841 row = rowidxs[i]; 842 col = colidxs[i]; 843 844 if( !SCIPsortedvecFindInt(&hescol[hesoffset[row]], col, hesoffset[row+1] - hesoffset[row], &pos) ) 845 { 846 SCIPerrorMessage("Could not find entry (%d, %d) in hessian sparsity\n", row, col); 847 return SCIP_ERROR; 848 } 849 850 values[hesoffset[row] + pos] += weight * h[i]; 851 } 852 853 return SCIP_OKAY; 854 } 855 856 /** prints a name, if available, makes sure it has not more than 64 characters, and adds a unique prefix if the longnames flag is set */ 857 static 858 void printName( 859 char* buffer, /**< buffer to print to, has to be not NULL and should be at least 65 bytes */ 860 char* name, /**< name, or NULL */ 861 int idx, /**< index of var or cons which the name corresponds to */ 862 char prefix, /**< a letter (typically 'x' or 'e') to distinguish variable and equation names, if names[idx] is not available */ 863 const char* suffix, /**< a suffer to add to the name, or NULL */ 864 SCIP_Bool longnames /**< whether prefixes for long names should be added */ 865 ) 866 { 867 assert(idx >= 0 && idx < 100000); /* to ensure that we do not exceed the size of the buffer */ 868 869 if( longnames ) 870 { 871 if( name != NULL ) 872 (void) SCIPsnprintf(buffer, 64, "%c%05d%.*s%s", prefix, idx, suffix != NULL ? (int)(57-strlen(suffix)) : 57, name, suffix ? suffix : ""); 873 else 874 (void) SCIPsnprintf(buffer, 64, "%c%05d", prefix, idx); 875 } 876 else 877 { 878 if( name != NULL ) 879 { 880 assert(strlen(name) + (suffix != NULL ? strlen(suffix) : 0) <= 64); 881 (void) SCIPsnprintf(buffer, 64, "%s%s", name, suffix != NULL ? suffix : ""); 882 } 883 else 884 { 885 assert(1 + 5 + (suffix != NULL ? strlen(suffix) : 0) <= 64); 886 (void) SCIPsnprintf(buffer, 64, "%c%d%s", prefix, idx, suffix != NULL ? suffix : ""); 887 } 888 } 889 } 890 891 /** prints a function */ 892 static 893 SCIP_RETCODE printFunction( 894 SCIP* scip, /**< SCIP data structure */ 895 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 896 FILE* file, /**< file to print to, has to be not NULL */ 897 SCIP_NLPIORACLECONS* cons, /**< constraint which function to print */ 898 SCIP_Bool longvarnames /**< whether variable names need to be shorten to 64 characters */ 899 ) 900 { /*lint --e{715}*/ 901 int i; 902 char namebuf[70]; 903 904 SCIPdebugMessage("%p print function\n", (void*)oracle); 905 906 assert(oracle != NULL); 907 assert(file != NULL); 908 assert(cons != NULL); 909 910 for( i = 0; i < cons->nlinidxs; ++i ) 911 { 912 printName(namebuf, oracle->varnames != NULL ? oracle->varnames[cons->linidxs[i]] : NULL, cons->linidxs[i], 'x', NULL, longvarnames); 913 SCIPinfoMessage(scip, file, "%+.15g*%s", cons->lincoefs[i], namebuf); 914 if( i % 10 == 9 ) 915 SCIPinfoMessage(scip, file, "\n"); 916 } 917 918 if( cons->expr != NULL ) 919 { 920 /* TODO SCIPprintExpr does not use the variable names in oracle->varnames, probably that should be changed */ 921 SCIPinfoMessage(scip, file, " +"); 922 SCIP_CALL( SCIPprintExpr(scip, cons->expr, file) ); 923 } 924 925 return SCIP_OKAY; 926 } 927 928 /** returns whether an expression contains nonsmooth operands (min, max, abs, ...) */ 929 static 930 SCIP_RETCODE exprIsNonSmooth( 931 SCIP* scip, /**< SCIP data structure */ 932 SCIP_EXPR* expr, /**< expression */ 933 SCIP_Bool* nonsmooth /**< buffer to store whether expression seems nonsmooth */ 934 ) 935 { 936 SCIP_EXPRITER* it; 937 938 assert(expr != NULL); 939 assert(nonsmooth != NULL); 940 941 *nonsmooth = FALSE; 942 943 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 944 SCIP_CALL( SCIPexpriterInit(it, expr, SCIP_EXPRITER_DFS, FALSE) ); 945 946 for( ; !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) 947 { 948 const char* hdlrname; 949 if( SCIPisExprSignpower(scip, expr) ) 950 { 951 *nonsmooth = TRUE; 952 break; 953 } 954 hdlrname = SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)); 955 if( strcmp(hdlrname, "abs") == 0 ) 956 { 957 *nonsmooth = TRUE; 958 break; 959 } 960 if( strcmp(hdlrname, "min") == 0 ) 961 { 962 *nonsmooth = TRUE; 963 break; 964 } 965 if( strcmp(hdlrname, "max") == 0 ) 966 { 967 *nonsmooth = TRUE; 968 break; 969 } 970 } 971 972 SCIPfreeExpriter(&it); 973 974 return SCIP_OKAY; 975 } 976 977 /**@} */ 978 979 /**@name public function */ 980 /**@{ */ 981 982 /** creates an NLPIORACLE data structure */ 983 SCIP_RETCODE SCIPnlpiOracleCreate( 984 SCIP* scip, /**< SCIP data structure */ 985 SCIP_NLPIORACLE** oracle /**< pointer to store NLPIORACLE data structure */ 986 ) 987 { 988 SCIP_Bool nlpieval; 989 990 assert(oracle != NULL); 991 992 SCIPdebugMessage("%p oracle create\n", (void*)oracle); 993 994 SCIP_CALL( SCIPallocMemory(scip, oracle) ); 995 BMSclearMemory(*oracle); 996 997 SCIPdebugMessage("Oracle initializes expression interpreter %s\n", SCIPexprintGetName()); 998 SCIP_CALL( SCIPexprintCreate(scip, &(*oracle)->exprinterpreter) ); 999 1000 SCIP_CALL( SCIPcreateClock(scip, &(*oracle)->evalclock) ); 1001 1002 SCIP_CALL( SCIPgetBoolParam(scip, "timing/nlpieval", &nlpieval) ); 1003 if( !nlpieval ) 1004 SCIPsetClockEnabled((*oracle)->evalclock, FALSE); 1005 1006 /* create zero objective function */ 1007 SCIP_CALL( createConstraint(scip, *oracle, &(*oracle)->objective, 0, NULL, NULL, NULL, 0.0, 0.0, NULL) ); 1008 1009 return SCIP_OKAY; 1010 } 1011 1012 /** frees an NLPIORACLE data structure */ 1013 SCIP_RETCODE SCIPnlpiOracleFree( 1014 SCIP* scip, /**< SCIP data structure */ 1015 SCIP_NLPIORACLE** oracle /**< pointer to NLPIORACLE data structure */ 1016 ) 1017 { 1018 assert(oracle != NULL); 1019 assert(*oracle != NULL); 1020 1021 SCIPdebugMessage("%p oracle free\n", (void*)oracle); 1022 1023 invalidateJacobiSparsity(scip, *oracle); 1024 invalidateHessianLagSparsity(scip, *oracle); 1025 1026 SCIP_CALL( freeConstraint(scip, *oracle, &(*oracle)->objective, FALSE) ); 1027 SCIP_CALL( freeConstraints(scip, *oracle) ); 1028 freeVariables(scip, *oracle); 1029 1030 SCIP_CALL( SCIPfreeClock(scip, &(*oracle)->evalclock) ); 1031 1032 SCIP_CALL( SCIPexprintFree(scip, &(*oracle)->exprinterpreter) ); 1033 1034 if( (*oracle)->name != NULL ) 1035 { 1036 SCIP_CALL( SCIPnlpiOracleSetProblemName(scip, *oracle, NULL) ); 1037 } 1038 1039 BMSfreeMemory(oracle); 1040 1041 return SCIP_OKAY; 1042 } 1043 1044 /** sets the problem name (used for printing) */ 1045 SCIP_RETCODE SCIPnlpiOracleSetProblemName( 1046 SCIP* scip, /**< SCIP data structure */ 1047 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1048 const char* name /**< name of problem */ 1049 ) 1050 { 1051 assert(oracle != NULL); 1052 1053 SCIPdebugMessage("%p set problem name\n", (void*)oracle); 1054 1055 if( oracle->name != NULL ) 1056 { 1057 SCIPfreeBlockMemoryArray(scip, &oracle->name, strlen(oracle->name)+1); 1058 } 1059 1060 if( name != NULL ) 1061 { 1062 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &oracle->name, name, strlen(name)+1) ); 1063 } 1064 1065 return SCIP_OKAY; 1066 } 1067 1068 /** gets the problem name, or NULL if none set */ 1069 const char* SCIPnlpiOracleGetProblemName( 1070 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1071 ) 1072 { 1073 assert(oracle != NULL); 1074 1075 SCIPdebugMessage("%p get problem name\n", (void*)oracle); 1076 1077 return oracle->name; 1078 } 1079 1080 /** adds variables */ 1081 SCIP_RETCODE SCIPnlpiOracleAddVars( 1082 SCIP* scip, /**< SCIP data structure */ 1083 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1084 int nvars, /**< number of variables to add */ 1085 const SCIP_Real* lbs, /**< array with lower bounds of new variables, or NULL if all -infinity */ 1086 const SCIP_Real* ubs, /**< array with upper bounds of new variables, or NULL if all +infinity */ 1087 const char** varnames /**< array with names of new variables, or NULL if no names should be stored */ 1088 ) 1089 { 1090 int i; 1091 1092 assert(oracle != NULL); 1093 1094 SCIPdebugMessage("%p add vars\n", (void*)oracle); 1095 1096 if( nvars == 0 ) 1097 return SCIP_OKAY; 1098 1099 assert(nvars > 0); 1100 1101 SCIP_CALL( ensureVarsSize(scip, oracle, oracle->nvars + nvars) ); 1102 1103 if( lbs != NULL ) 1104 { 1105 BMScopyMemoryArray(&oracle->varlbs[oracle->nvars], lbs, nvars); 1106 } 1107 else 1108 for( i = 0; i < nvars; ++i ) 1109 oracle->varlbs[oracle->nvars+i] = -SCIPinfinity(scip); 1110 1111 if( ubs != NULL ) 1112 { 1113 BMScopyMemoryArray(&oracle->varubs[oracle->nvars], ubs, nvars); 1114 1115 /* ensure variable bounds are consistent */ 1116 for( i = oracle->nvars; i < oracle->nvars + nvars; ++i ) 1117 { 1118 if( oracle->varlbs[i] > oracle->varubs[i] ) 1119 { 1120 assert(EPSEQ(oracle->varlbs[i], oracle->varubs[i], SCIP_DEFAULT_EPSILON)); 1121 oracle->varlbs[i] = oracle->varubs[i]; 1122 } 1123 } 1124 } 1125 else 1126 for( i = 0; i < nvars; ++i ) 1127 oracle->varubs[oracle->nvars+i] = SCIPinfinity(scip); 1128 1129 if( varnames != NULL ) 1130 { 1131 if( oracle->varnames == NULL ) 1132 { 1133 SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &oracle->varnames, oracle->varssize) ); 1134 } 1135 1136 for( i = 0; i < nvars; ++i ) 1137 { 1138 if( varnames[i] != NULL ) 1139 { 1140 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &oracle->varnames[oracle->nvars+i], varnames[i], strlen(varnames[i])+1) ); 1141 } 1142 else 1143 oracle->varnames[oracle->nvars+i] = NULL; 1144 } 1145 } 1146 else if( oracle->varnames != NULL ) 1147 { 1148 BMSclearMemoryArray(&oracle->varnames[oracle->nvars], nvars); 1149 } 1150 1151 BMSclearMemoryArray(&oracle->varlincount[oracle->nvars], nvars); 1152 BMSclearMemoryArray(&oracle->varnlcount[oracle->nvars], nvars); 1153 1154 /* @TODO update sparsity pattern by extending heslagoffsets */ 1155 invalidateHessianLagSparsity(scip, oracle); 1156 1157 oracle->nvars += nvars; 1158 1159 return SCIP_OKAY; 1160 } 1161 1162 /** adds constraints 1163 * 1164 * linear coefficients: row(=constraint) oriented matrix; 1165 * quadratic coefficients: row oriented matrix for each constraint 1166 */ 1167 SCIP_RETCODE SCIPnlpiOracleAddConstraints( 1168 SCIP* scip, /**< SCIP data structure */ 1169 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1170 int nconss, /**< number of constraints to add */ 1171 const SCIP_Real* lhss, /**< array with left-hand sides of constraints, or NULL if all -infinity */ 1172 const SCIP_Real* rhss, /**< array with right-hand sides of constraints, or NULL if all +infinity */ 1173 const int* nlininds, /**< number of linear coefficients for each constraint, may be NULL in case of no linear part */ 1174 int* const* lininds, /**< indices of variables for linear coefficients for each constraint, may be NULL in case of no linear part */ 1175 SCIP_Real* const* linvals, /**< values of linear coefficient for each constraint, may be NULL in case of no linear part */ 1176 SCIP_EXPR** exprs, /**< NULL if no nonlinear parts, otherwise exprs[.] gives nonlinear part, 1177 * or NULL if no nonlinear part in this constraint */ 1178 const char** consnames /**< names of new constraints, or NULL if no names should be stored */ 1179 ) 1180 { /*lint --e{715}*/ 1181 SCIP_NLPIORACLECONS* cons; 1182 SCIP_Bool addednlcon; /* whether a nonlinear constraint was added */ 1183 int c; 1184 1185 assert(oracle != NULL); 1186 1187 SCIPdebugMessage("%p add constraints\n", (void*)oracle); 1188 1189 if( nconss == 0 ) 1190 return SCIP_OKAY; 1191 1192 assert(nconss > 0); 1193 1194 addednlcon = FALSE; 1195 1196 invalidateJacobiSparsity(scip, oracle); /* @TODO we could also update (extend) the sparsity pattern */ 1197 1198 SCIP_CALL( ensureConssSize(scip, oracle, oracle->nconss + nconss) ); 1199 for( c = 0; c < nconss; ++c ) 1200 { 1201 SCIP_CALL( createConstraint(scip, oracle, &cons, 1202 nlininds != NULL ? nlininds[c] : 0, 1203 lininds != NULL ? lininds[c] : NULL, 1204 linvals != NULL ? linvals[c] : NULL, 1205 exprs != NULL ? exprs[c] : NULL, 1206 lhss != NULL ? lhss[c] : -SCIPinfinity(scip), 1207 rhss != NULL ? rhss[c] : SCIPinfinity(scip), 1208 consnames != NULL ? consnames[c] : NULL 1209 ) ); 1210 1211 if( cons->expr != NULL ) 1212 addednlcon = TRUE; 1213 1214 oracle->conss[oracle->nconss+c] = cons; 1215 } 1216 oracle->nconss += nconss; 1217 1218 if( addednlcon == TRUE ) 1219 invalidateHessianLagSparsity(scip, oracle); 1220 1221 return SCIP_OKAY; 1222 } 1223 1224 /** sets or overwrites objective, a minimization problem is expected 1225 * 1226 * May change sparsity pattern. 1227 */ 1228 SCIP_RETCODE SCIPnlpiOracleSetObjective( 1229 SCIP* scip, /**< SCIP data structure */ 1230 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1231 const SCIP_Real constant, /**< constant part of objective */ 1232 int nlin, /**< number of linear variable coefficients */ 1233 const int* lininds, /**< indices of linear variables, or NULL if no linear part */ 1234 const SCIP_Real* linvals, /**< coefficients of linear variables, or NULL if no linear part */ 1235 SCIP_EXPR* expr /**< expression of nonlinear part, or NULL if no nonlinear part */ 1236 ) 1237 { /*lint --e{715}*/ 1238 assert(oracle != NULL); 1239 assert(!SCIPisInfinity(scip, REALABS(constant))); 1240 1241 SCIPdebugMessage("%p set objective\n", (void*)oracle); 1242 1243 if( expr != NULL || oracle->objective->expr != NULL ) 1244 invalidateHessianLagSparsity(scip, oracle); 1245 1246 /* clear previous objective */ 1247 SCIP_CALL( freeConstraint(scip, oracle, &oracle->objective, TRUE) ); 1248 1249 /* create new objective */ 1250 SCIP_CALL( createConstraint(scip, oracle, &oracle->objective, 1251 nlin, lininds, linvals, expr, constant, constant, NULL) ); 1252 1253 return SCIP_OKAY; 1254 } 1255 1256 /** change variable bounds */ 1257 SCIP_RETCODE SCIPnlpiOracleChgVarBounds( 1258 SCIP* scip, /**< SCIP data structure */ 1259 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1260 int nvars, /**< number of variables to change bounds */ 1261 const int* indices, /**< indices of variables to change bounds */ 1262 const SCIP_Real* lbs, /**< new lower bounds, or NULL if all should be -infty */ 1263 const SCIP_Real* ubs /**< new upper bounds, or NULL if all should be +infty */ 1264 ) 1265 { 1266 int i; 1267 1268 assert(oracle != NULL); 1269 assert(indices != NULL || nvars == 0); 1270 1271 SCIPdebugMessage("%p chg var bounds\n", (void*)oracle); 1272 1273 for( i = 0; i < nvars; ++i ) 1274 { 1275 assert(indices != NULL); 1276 assert(indices[i] >= 0); 1277 assert(indices[i] < oracle->nvars); 1278 1279 oracle->varlbs[indices[i]] = (lbs != NULL ? lbs[i] : -SCIPinfinity(scip)); 1280 oracle->varubs[indices[i]] = (ubs != NULL ? ubs[i] : SCIPinfinity(scip)); 1281 1282 if( oracle->varlbs[indices[i]] > oracle->varubs[indices[i]] ) 1283 { 1284 /* inconsistent bounds; let's assume it's due to rounding and make them equal */ 1285 assert(EPSEQ(oracle->varlbs[indices[i]], oracle->varubs[indices[i]], SCIP_DEFAULT_EPSILON)); 1286 oracle->varlbs[indices[i]] = oracle->varubs[indices[i]]; 1287 } 1288 } 1289 1290 return SCIP_OKAY; 1291 } 1292 1293 /** change constraint sides */ 1294 SCIP_RETCODE SCIPnlpiOracleChgConsSides( 1295 SCIP* scip, /**< SCIP data structure */ 1296 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1297 int nconss, /**< number of constraints to change bounds */ 1298 const int* indices, /**< indices of constraints to change bounds */ 1299 const SCIP_Real* lhss, /**< new left-hand sides, or NULL if all should be -infty */ 1300 const SCIP_Real* rhss /**< new right-hand sides, or NULL if all should be +infty */ 1301 ) 1302 { 1303 int i; 1304 1305 assert(oracle != NULL); 1306 assert(indices != NULL || nconss == 0); 1307 1308 SCIPdebugMessage("%p chg cons sides\n", (void*)oracle); 1309 1310 for( i = 0; i < nconss; ++i ) 1311 { 1312 assert(indices != NULL); 1313 assert(indices[i] >= 0); 1314 assert(indices[i] < oracle->nconss); 1315 1316 oracle->conss[indices[i]]->lhs = (lhss != NULL ? lhss[i] : -SCIPinfinity(scip)); 1317 oracle->conss[indices[i]]->rhs = (rhss != NULL ? rhss[i] : SCIPinfinity(scip)); 1318 if( oracle->conss[indices[i]]->lhs > oracle->conss[indices[i]]->rhs ) 1319 { 1320 assert(EPSEQ(oracle->conss[indices[i]]->lhs, oracle->conss[indices[i]]->rhs, SCIP_DEFAULT_EPSILON)); 1321 oracle->conss[indices[i]]->lhs = oracle->conss[indices[i]]->rhs; 1322 } 1323 } 1324 1325 return SCIP_OKAY; 1326 } 1327 1328 /** deletes a set of variables */ 1329 SCIP_RETCODE SCIPnlpiOracleDelVarSet( 1330 SCIP* scip, /**< SCIP data structure */ 1331 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1332 int* delstats /**< deletion status of vars in input (1 if var should be deleted, 0 if not); 1333 * new position of var in output (-1 if var was deleted) */ 1334 ) 1335 { /*lint --e{715}*/ 1336 int c; 1337 int lastgood; /* index of the last variable that should be kept */ 1338 SCIP_NLPIORACLECONS* cons; 1339 SCIP_EXPRITER* it; 1340 1341 assert(oracle != NULL); 1342 1343 SCIPdebugMessage("%p del var set\n", (void*)oracle); 1344 1345 invalidateJacobiSparsity(scip, oracle); 1346 invalidateHessianLagSparsity(scip, oracle); 1347 1348 lastgood = oracle->nvars - 1; 1349 while( lastgood >= 0 && delstats[lastgood] == 1 ) 1350 --lastgood; 1351 if( lastgood < 0 ) 1352 { 1353 /* all variables should be deleted */ 1354 assert(oracle->nconss == 0); /* we could relax this by checking that all constraints are constant */ 1355 oracle->objective->nlinidxs = 0; 1356 for( c = 0; c < oracle->nvars; ++c ) 1357 delstats[c] = -1; 1358 freeVariables(scip, oracle); 1359 return SCIP_OKAY; 1360 } 1361 1362 /* delete variables at the end */ 1363 for( c = oracle->nvars - 1; c > lastgood; --c ) 1364 { 1365 if( oracle->varnames && oracle->varnames[c] != NULL ) 1366 { 1367 SCIPfreeBlockMemoryArray(scip, &oracle->varnames[c], strlen(oracle->varnames[c])+1); 1368 } 1369 delstats[c] = -1; 1370 } 1371 1372 /* go through variables from the beginning on 1373 * if variable should be deleted, free it and move lastgood variable to this position 1374 * then update lastgood */ 1375 for( c = 0; c <= lastgood; ++c ) 1376 { 1377 if( delstats[c] == 0 ) 1378 { /* variable should not be deleted and is kept on position c */ 1379 delstats[c] = c; 1380 continue; 1381 } 1382 assert(delstats[c] == 1); /* variable should be deleted */ 1383 1384 if( oracle->varnames != NULL && oracle->varnames[c] != NULL ) 1385 { 1386 SCIPfreeBlockMemoryArray(scip, &oracle->varnames[c], strlen(oracle->varnames[c])+1); 1387 } 1388 delstats[c] = -1; 1389 1390 /* move variable at position lastgood to position c */ 1391 SCIP_CALL( moveVariable(scip, oracle, lastgood, c) ); 1392 delstats[lastgood] = c; /* mark that lastgood variable is now at position c */ 1393 1394 /* move lastgood forward, delete variables on the way */ 1395 --lastgood; 1396 while( lastgood > c && delstats[lastgood] == 1) 1397 { 1398 if( oracle->varnames && oracle->varnames[lastgood] != NULL ) 1399 { 1400 SCIPfreeBlockMemoryArray(scip, &oracle->varnames[lastgood], strlen(oracle->varnames[lastgood])+1); 1401 } 1402 delstats[lastgood] = -1; 1403 --lastgood; 1404 } 1405 } 1406 assert(c == lastgood); 1407 1408 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 1409 1410 for( c = -1; c < oracle->nconss; ++c ) 1411 { 1412 cons = c < 0 ? oracle->objective : oracle->conss[c]; 1413 assert(cons != NULL); 1414 1415 /* update indices in linear part, sort indices, and then clear elements that are marked as deleted */ 1416 mapIndices(delstats, cons->nlinidxs, cons->linidxs); 1417 SCIPsortIntReal(cons->linidxs, cons->lincoefs, cons->nlinidxs); 1418 clearDeletedLinearElements(&cons->linidxs, &cons->lincoefs, &cons->nlinidxs); 1419 1420 if( cons->expr != NULL ) 1421 { 1422 /* update variable indices in varidx expressions */ 1423 SCIP_EXPR* expr; 1424 SCIP_Bool keptvar = FALSE; /* whether any of the variables in expr was not deleted */ 1425 #ifndef NDEBUG 1426 SCIP_Bool delvar = FALSE; /* whether any of the variables in expr was deleted */ 1427 #endif 1428 1429 SCIP_CALL( SCIPexpriterInit(it, cons->expr, SCIP_EXPRITER_DFS, FALSE) ); 1430 for( expr = cons->expr; !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) 1431 { 1432 if( !SCIPisExprVaridx(scip, expr) ) 1433 continue; 1434 1435 if( delstats[SCIPgetIndexExprVaridx(expr)] >= 0 ) 1436 { 1437 /* if variable is not deleted, then set its new index */ 1438 keptvar = TRUE; 1439 SCIPsetIndexExprVaridx(expr, delstats[SCIPgetIndexExprVaridx(expr)]); 1440 1441 /* if variable is kept, then there must not have been any variable that was deleted */ 1442 assert(!delvar); 1443 } 1444 else 1445 { 1446 #ifndef NDEBUG 1447 delvar = TRUE; 1448 #endif 1449 /* if variable is deleted, then there must not have been any variable that was kept 1450 * (either all variables are deleted, which removes the expr, or none) 1451 */ 1452 assert(!keptvar); 1453 } 1454 } 1455 if( !keptvar ) 1456 { 1457 SCIP_CALL( SCIPexprintFreeData(scip, oracle->exprinterpreter, cons->expr, &cons->exprintdata) ); 1458 SCIP_CALL( SCIPreleaseExpr(scip, &cons->expr) ); 1459 } 1460 } 1461 } 1462 1463 SCIPfreeExpriter(&it); 1464 1465 oracle->nvars = lastgood+1; 1466 1467 return SCIP_OKAY; 1468 } 1469 1470 /** deletes a set of constraints */ 1471 SCIP_RETCODE SCIPnlpiOracleDelConsSet( 1472 SCIP* scip, /**< SCIP data structure */ 1473 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1474 int* delstats /**< array with deletion status of rows in input (1 if row should be deleted, 0 if not); 1475 * new position of row in output (-1 if row was deleted) */ 1476 ) 1477 { /*lint --e{715}*/ 1478 int c; 1479 int lastgood; /* index of the last constraint that should be kept */ 1480 1481 assert(oracle != NULL); 1482 1483 SCIPdebugMessage("%p del cons set\n", (void*)oracle); 1484 1485 invalidateJacobiSparsity(scip, oracle); 1486 invalidateHessianLagSparsity(scip, oracle); 1487 1488 lastgood = oracle->nconss - 1; 1489 while( lastgood >= 0 && delstats[lastgood] == 1) 1490 --lastgood; 1491 if( lastgood < 0 ) 1492 { 1493 /* all constraints should be deleted */ 1494 for( c = 0; c < oracle->nconss; ++c ) 1495 delstats[c] = -1; 1496 SCIP_CALL( freeConstraints(scip, oracle) ); 1497 1498 /* the previous call did not keep variable counts uptodate 1499 * since we only have an objective function left, we reset the counts to the ones of the objective 1500 */ 1501 BMSclearMemoryArray(oracle->varlincount, oracle->nvars); 1502 BMSclearMemoryArray(oracle->varnlcount, oracle->nvars); 1503 SCIP_CALL( updateVariableCounts(scip, oracle, 1, oracle->objective->nlinidxs, oracle->objective->linidxs, oracle->objective->expr) ); 1504 1505 return SCIP_OKAY; 1506 } 1507 1508 /* delete constraints at the end */ 1509 for( c = oracle->nconss - 1; c > lastgood; --c ) 1510 { 1511 SCIP_CALL( freeConstraint(scip, oracle, &oracle->conss[c], TRUE) ); 1512 assert(oracle->conss[c] == NULL); 1513 delstats[c] = -1; 1514 } 1515 1516 /* go through constraint from the beginning on 1517 * if constraint should be deleted, free it and move lastgood constraint to this position 1518 * then update lastgood */ 1519 for( c = 0; c <= lastgood; ++c ) 1520 { 1521 if( delstats[c] == 0 ) 1522 { 1523 /* constraint should not be deleted and is kept on position c */ 1524 delstats[c] = c; 1525 continue; 1526 } 1527 assert(delstats[c] == 1); /* constraint should be deleted */ 1528 1529 SCIP_CALL( freeConstraint(scip, oracle, &oracle->conss[c], TRUE) ); 1530 assert(oracle->conss[c] == NULL); 1531 delstats[c] = -1; 1532 1533 /* move constraint at position lastgood to position c */ 1534 oracle->conss[c] = oracle->conss[lastgood]; 1535 assert(oracle->conss[c] != NULL); 1536 delstats[lastgood] = c; /* mark that lastgood constraint is now at position c */ 1537 oracle->conss[lastgood] = NULL; 1538 --lastgood; 1539 1540 /* move lastgood forward, delete constraints on the way */ 1541 while( lastgood > c && delstats[lastgood] == 1) 1542 { 1543 SCIP_CALL( freeConstraint(scip, oracle, &oracle->conss[lastgood], TRUE) ); 1544 assert(oracle->conss[lastgood] == NULL); 1545 delstats[lastgood] = -1; 1546 --lastgood; 1547 } 1548 } 1549 assert(c == lastgood+1); 1550 1551 oracle->nconss = lastgood+1; 1552 1553 return SCIP_OKAY; 1554 } 1555 1556 /** changes (or adds) linear coefficients in one constraint or objective */ 1557 SCIP_RETCODE SCIPnlpiOracleChgLinearCoefs( 1558 SCIP* scip, /**< SCIP data structure */ 1559 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1560 int considx, /**< index of constraint where linear coefficients should be changed, or -1 for objective */ 1561 int nentries, /**< number of coefficients to change */ 1562 const int* varidxs, /**< array with indices of variables which coefficients should be changed */ 1563 const SCIP_Real* newcoefs /**< array with new coefficients of variables */ 1564 ) 1565 { /*lint --e{715}*/ 1566 SCIP_NLPIORACLECONS* cons; 1567 SCIP_Bool needsort; 1568 int i; 1569 1570 SCIPdebugMessage("%p chg linear coefs\n", (void*)oracle); 1571 1572 assert(oracle != NULL); 1573 assert(varidxs != NULL || nentries == 0); 1574 assert(newcoefs != NULL || nentries == 0); 1575 assert(considx >= -1); 1576 assert(considx < oracle->nconss); 1577 1578 if( nentries == 0 ) 1579 return SCIP_OKAY; 1580 1581 SCIPdebugMessage("change %d linear coefficients in cons %d\n", nentries, considx); 1582 1583 needsort = FALSE; 1584 1585 cons = considx < 0 ? oracle->objective : oracle->conss[considx]; 1586 1587 if( cons->linsize == 0 ) 1588 { 1589 /* first time we have linear coefficients in this constraint (or objective) */ 1590 assert(cons->linidxs == NULL); 1591 assert(cons->lincoefs == NULL); 1592 1593 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &cons->linidxs, varidxs, nentries) ); 1594 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &cons->lincoefs, newcoefs, nentries) ); 1595 cons->linsize = nentries; 1596 cons->nlinidxs = nentries; 1597 1598 SCIP_CALL( updateVariableCounts(scip, oracle, 1, nentries, varidxs, NULL) ); 1599 1600 needsort = TRUE; 1601 } 1602 else 1603 { 1604 int pos; 1605 1606 for( i = 0; i < nentries; ++i ) 1607 { 1608 assert(varidxs[i] >= 0); /*lint !e613*/ 1609 assert(varidxs[i] < oracle->nvars); /*lint !e613*/ 1610 1611 if( SCIPsortedvecFindInt(cons->linidxs, varidxs[i], cons->nlinidxs, &pos) ) /*lint !e613*/ 1612 { 1613 SCIPdebugMessage("replace coefficient of var %d at pos %d by %g\n", varidxs[i], pos, newcoefs[i]); /*lint !e613*/ 1614 1615 cons->lincoefs[pos] = newcoefs[i]; /*lint !e613*/ 1616 1617 /* remember that we need to sort/merge/squeeze array if coefficient became zero here */ 1618 needsort |= (newcoefs[i] == 0.0); /*lint !e613 !e514*/ 1619 1620 if( newcoefs[i] == 0.0 ) 1621 { 1622 --oracle->varlincount[varidxs[i]]; 1623 assert(oracle->varlincount[varidxs[i]] >= 0); 1624 } 1625 } 1626 else if( newcoefs[i] != 0.0 ) /*lint !e613*/ 1627 { 1628 /* append new entry */ 1629 SCIPdebugMessage("add coefficient of var %d at pos %d, value %g\n", varidxs[i], cons->nlinidxs, newcoefs[i]); /*lint !e613*/ 1630 1631 SCIP_CALL( ensureConsLinSize(scip, cons, cons->nlinidxs + (nentries-i)) ); 1632 cons->linidxs[cons->nlinidxs] = varidxs[i]; /*lint !e613*/ 1633 cons->lincoefs[cons->nlinidxs] = newcoefs[i]; /*lint !e613*/ 1634 ++cons->nlinidxs; 1635 1636 ++oracle->varlincount[varidxs[i]]; 1637 1638 needsort = TRUE; 1639 } 1640 } 1641 } 1642 1643 if( needsort ) 1644 { 1645 invalidateJacobiSparsity(scip, oracle); 1646 sortLinearCoefficients(&cons->nlinidxs, cons->linidxs, cons->lincoefs); 1647 } 1648 1649 return SCIP_OKAY; 1650 } 1651 1652 /** replaces expression of one constraint or objective */ 1653 SCIP_RETCODE SCIPnlpiOracleChgExpr( 1654 SCIP* scip, /**< SCIP data structure */ 1655 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1656 int considx, /**< index of constraint where expression should be changed, or -1 for objective */ 1657 SCIP_EXPR* expr /**< new expression, or NULL */ 1658 ) 1659 { 1660 SCIP_NLPIORACLECONS* cons; 1661 1662 SCIPdebugMessage("%p chg expr\n", (void*)oracle); 1663 1664 assert(oracle != NULL); 1665 assert(considx >= -1); 1666 assert(considx < oracle->nconss); 1667 1668 invalidateHessianLagSparsity(scip, oracle); 1669 invalidateJacobiSparsity(scip, oracle); 1670 1671 cons = considx < 0 ? oracle->objective : oracle->conss[considx]; 1672 1673 /* free previous expression */ 1674 if( cons->expr != NULL ) 1675 { 1676 SCIP_CALL( updateVariableCounts(scip, oracle, -1, 0, NULL, cons->expr) ); 1677 SCIP_CALL( SCIPexprintFreeData(scip, oracle->exprinterpreter, cons->expr, &cons->exprintdata) ); 1678 SCIP_CALL( SCIPreleaseExpr(scip, &cons->expr) ); 1679 } 1680 1681 /* if user did not want to set new expr, then we are done */ 1682 if( expr == NULL ) 1683 return SCIP_OKAY; 1684 1685 assert(oracle->exprinterpreter != NULL); 1686 1687 /* install new expression */ 1688 cons->expr = expr; 1689 SCIPcaptureExpr(cons->expr); 1690 SCIP_CALL( SCIPexprintCompile(scip, oracle->exprinterpreter, cons->expr, &cons->exprintdata) ); 1691 1692 /* keep variable counts up to date */ 1693 SCIP_CALL( updateVariableCounts(scip, oracle, 1, 0, NULL, cons->expr) ); 1694 1695 return SCIP_OKAY; 1696 } 1697 1698 /** changes the constant value in the objective function */ /*lint -e{715}*/ 1699 SCIP_RETCODE SCIPnlpiOracleChgObjConstant( 1700 SCIP* scip, /**< SCIP data structure */ 1701 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1702 SCIP_Real objconstant /**< new value for objective constant */ 1703 ) 1704 { /*lint --e{715}*/ 1705 assert(oracle != NULL); 1706 1707 SCIPdebugMessage("%p chg obj constant\n", (void*)oracle); 1708 1709 oracle->objective->lhs = objconstant; 1710 oracle->objective->rhs = objconstant; 1711 1712 return SCIP_OKAY; 1713 } 1714 1715 /** gives the current number of variables */ 1716 int SCIPnlpiOracleGetNVars( 1717 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1718 ) 1719 { 1720 assert(oracle != NULL); 1721 1722 return oracle->nvars; 1723 } 1724 1725 /** gives the current number of constraints */ 1726 int SCIPnlpiOracleGetNConstraints( 1727 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1728 ) 1729 { 1730 assert(oracle != NULL); 1731 1732 return oracle->nconss; 1733 } 1734 1735 /** gives the variables lower bounds */ 1736 const SCIP_Real* SCIPnlpiOracleGetVarLbs( 1737 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1738 ) 1739 { 1740 assert(oracle != NULL); 1741 1742 return oracle->varlbs; 1743 } 1744 1745 /** gives the variables upper bounds */ 1746 const SCIP_Real* SCIPnlpiOracleGetVarUbs( 1747 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1748 ) 1749 { 1750 assert(oracle != NULL); 1751 1752 return oracle->varubs; 1753 } 1754 1755 /** gives the variables names, or NULL if not set */ 1756 char** SCIPnlpiOracleGetVarNames( 1757 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1758 ) 1759 { 1760 assert(oracle != NULL); 1761 1762 return oracle->varnames; 1763 } 1764 1765 /** indicates whether variable appears nonlinear in any objective or constraint */ /*lint --e{715}*/ 1766 SCIP_Bool SCIPnlpiOracleIsVarNonlinear( 1767 SCIP* scip, /**< SCIP data structure */ 1768 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1769 int varidx /**< the variable to check */ 1770 ) 1771 { 1772 assert(oracle != NULL); 1773 assert(varidx >= 0); 1774 assert(varidx < oracle->nvars); 1775 assert(oracle->varnlcount != NULL); 1776 1777 return oracle->varnlcount[varidx] > 0; 1778 } 1779 1780 /** returns number of linear and nonlinear appearances of variables in objective and constraints */ /*lint --e{715}*/ 1781 void SCIPnlpiOracleGetVarCounts( 1782 SCIP* scip, /**< SCIP data structure */ 1783 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1784 const int** lincounts, /**< buffer to return pointer to array of counts of linear appearances */ 1785 const int** nlcounts /**< buffer to return pointer to array of counts of nonlinear appearances */ 1786 ) 1787 { 1788 assert(oracle != NULL); 1789 assert(lincounts != NULL); 1790 assert(nlcounts != NULL); 1791 1792 *lincounts = oracle->varlincount; 1793 *nlcounts = oracle->varnlcount; 1794 } 1795 1796 /** gives constant term of objective */ 1797 SCIP_Real SCIPnlpiOracleGetObjectiveConstant( 1798 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1799 ) 1800 { 1801 assert(oracle != NULL); 1802 assert(oracle->objective->lhs == oracle->objective->rhs); /*lint !e777*/ 1803 1804 return oracle->objective->lhs; 1805 } 1806 1807 /** gives left-hand side of a constraint */ 1808 SCIP_Real SCIPnlpiOracleGetConstraintLhs( 1809 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1810 int considx /**< constraint index */ 1811 ) 1812 { 1813 assert(oracle != NULL); 1814 assert(considx >= 0); 1815 assert(considx < oracle->nconss); 1816 1817 return oracle->conss[considx]->lhs; 1818 } 1819 1820 /** gives right-hand side of a constraint */ 1821 SCIP_Real SCIPnlpiOracleGetConstraintRhs( 1822 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1823 int considx /**< constraint index */ 1824 ) 1825 { 1826 assert(oracle != NULL); 1827 assert(considx >= 0); 1828 assert(considx < oracle->nconss); 1829 1830 return oracle->conss[considx]->rhs; 1831 } 1832 1833 /** gives name of a constraint, may be NULL */ 1834 char* SCIPnlpiOracleGetConstraintName( 1835 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1836 int considx /**< constraint index */ 1837 ) 1838 { 1839 assert(oracle != NULL); 1840 assert(considx >= 0); 1841 assert(considx < oracle->nconss); 1842 1843 return oracle->conss[considx]->name; 1844 } 1845 1846 /** indicates whether constraint is nonlinear */ 1847 SCIP_Bool SCIPnlpiOracleIsConstraintNonlinear( 1848 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1849 int considx /**< index of constraint for which nonlinearity status is returned, or -1 for objective */ 1850 ) 1851 { 1852 SCIP_NLPIORACLECONS* cons; 1853 1854 assert(oracle != NULL); 1855 assert(considx >= -1); 1856 assert(considx < oracle->nconss); 1857 1858 cons = considx < 0 ? oracle->objective : oracle->conss[considx]; 1859 1860 return cons->expr != NULL; 1861 } 1862 1863 /** gives the evaluation capabilities that are shared among all expressions in the problem */ 1864 SCIP_EXPRINTCAPABILITY SCIPnlpiOracleGetEvalCapability( 1865 SCIP* scip, /**< SCIP data structure */ 1866 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 1867 ) 1868 { 1869 int c; 1870 SCIP_EXPRINTCAPABILITY evalcapability; 1871 1872 assert(oracle != NULL); 1873 1874 if( oracle->objective->expr != NULL ) 1875 evalcapability = SCIPexprintGetExprCapability(scip, oracle->exprinterpreter, oracle->objective->expr, oracle->objective->exprintdata); 1876 else 1877 evalcapability = SCIP_EXPRINTCAPABILITY_ALL; 1878 1879 for( c = 0; c < oracle->nconss; ++c ) 1880 if( oracle->conss[c]->expr != NULL ) 1881 evalcapability &= SCIPexprintGetExprCapability(scip, oracle->exprinterpreter, oracle->conss[c]->expr, oracle->conss[c]->exprintdata); 1882 1883 return evalcapability; 1884 } 1885 1886 /** evaluates the objective function in a given point */ 1887 SCIP_RETCODE SCIPnlpiOracleEvalObjectiveValue( 1888 SCIP* scip, /**< SCIP data structure */ 1889 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1890 const SCIP_Real* x, /**< point where to evaluate */ 1891 SCIP_Real* objval /**< pointer to store objective value */ 1892 ) 1893 { 1894 SCIP_RETCODE retcode; 1895 1896 assert(oracle != NULL); 1897 1898 SCIPdebugMessage("%p eval obj value\n", (void*)oracle); 1899 1900 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 1901 retcode = evalFunctionValue(scip, oracle, oracle->objective, x, objval); 1902 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 1903 1904 assert(oracle->objective->lhs == oracle->objective->rhs); /*lint !e777*/ 1905 if( retcode == SCIP_OKAY ) 1906 *objval += oracle->objective->lhs; 1907 1908 return retcode; 1909 } 1910 1911 /** evaluates one constraint function in a given point */ 1912 SCIP_RETCODE SCIPnlpiOracleEvalConstraintValue( 1913 SCIP* scip, /**< SCIP data structure */ 1914 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1915 int considx, /**< index of constraint to evaluate */ 1916 const SCIP_Real* x, /**< point where to evaluate */ 1917 SCIP_Real* conval /**< pointer to store constraint value */ 1918 ) 1919 { 1920 SCIP_RETCODE retcode; 1921 1922 assert(oracle != NULL); 1923 assert(x != NULL || oracle->nvars == 0); 1924 assert(conval != NULL); 1925 1926 SCIPdebugMessage("%p eval cons value\n", (void*)oracle); 1927 1928 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 1929 retcode = evalFunctionValue(scip, oracle, oracle->conss[considx], x, conval); 1930 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 1931 1932 return retcode; 1933 } 1934 1935 /** evaluates all constraint functions in a given point */ 1936 SCIP_RETCODE SCIPnlpiOracleEvalConstraintValues( 1937 SCIP* scip, /**< SCIP data structure */ 1938 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1939 const SCIP_Real* x, /**< point where to evaluate */ 1940 SCIP_Real* convals /**< buffer to store constraint values */ 1941 ) 1942 { 1943 SCIP_RETCODE retcode = SCIP_OKAY; 1944 int i; 1945 1946 SCIPdebugMessage("%p eval cons values\n", (void*)oracle); 1947 1948 assert(oracle != NULL); 1949 assert(x != NULL || oracle->nvars == 0); 1950 assert(convals != NULL); 1951 1952 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 1953 for( i = 0; i < oracle->nconss; ++i ) 1954 { 1955 retcode = evalFunctionValue(scip, oracle, oracle->conss[i], x, &convals[i]); 1956 if( retcode != SCIP_OKAY ) 1957 break; 1958 } 1959 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 1960 1961 return retcode; 1962 } 1963 1964 /** computes the objective gradient in a given point 1965 * 1966 * @return SCIP_INVALIDDATA, if the function or its gradient could not be evaluated (domain error, etc.) 1967 */ 1968 SCIP_RETCODE SCIPnlpiOracleEvalObjectiveGradient( 1969 SCIP* scip, /**< SCIP data structure */ 1970 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 1971 const SCIP_Real* x, /**< point where to evaluate */ 1972 SCIP_Bool isnewx, /**< has the point x changed since the last call to some evaluation function? */ 1973 SCIP_Real* objval, /**< pointer to store objective value */ 1974 SCIP_Real* objgrad /**< pointer to store (dense) objective gradient */ 1975 ) 1976 { 1977 SCIP_RETCODE retcode; 1978 assert(oracle != NULL); 1979 1980 SCIPdebugMessage("%p eval obj grad\n", (void*)oracle); 1981 1982 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 1983 retcode = evalFunctionGradient(scip, oracle, oracle->objective, x, isnewx, objval, objgrad); 1984 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 1985 1986 assert(oracle->objective->lhs == oracle->objective->rhs); /*lint !e777*/ 1987 if( retcode == SCIP_OKAY ) 1988 *objval += oracle->objective->lhs; 1989 1990 return retcode; 1991 } 1992 1993 /** computes a constraints gradient in a given point 1994 * 1995 * @return SCIP_INVALIDDATA, if the function or its gradient could not be evaluated (domain error, etc.) 1996 */ 1997 SCIP_RETCODE SCIPnlpiOracleEvalConstraintGradient( 1998 SCIP* scip, /**< SCIP data structure */ 1999 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 2000 const int considx, /**< index of constraint to compute gradient for */ 2001 const SCIP_Real* x, /**< point where to evaluate */ 2002 SCIP_Bool isnewx, /**< has the point x changed since the last call to some evaluation function? */ 2003 SCIP_Real* conval, /**< pointer to store constraint value */ 2004 SCIP_Real* congrad /**< pointer to store (dense) constraint gradient */ 2005 ) 2006 { 2007 SCIP_RETCODE retcode; 2008 2009 assert(oracle != NULL); 2010 assert(x != NULL || oracle->nvars == 0); 2011 assert(conval != NULL); 2012 2013 SCIPdebugMessage("%p eval cons grad\n", (void*)oracle); 2014 2015 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 2016 retcode = evalFunctionGradient(scip, oracle, oracle->conss[considx], x, isnewx, conval, congrad); 2017 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 2018 2019 return retcode; 2020 } 2021 2022 /** gets sparsity pattern (rowwise) of Jacobian matrix 2023 * 2024 * Note that internal data is returned in *offset and *col, thus the user does not need to allocate memory there. 2025 * Adding or deleting constraints destroys the sparsity structure and make another call to this function necessary. 2026 */ 2027 SCIP_RETCODE SCIPnlpiOracleGetJacobianSparsity( 2028 SCIP* scip, /**< SCIP data structure */ 2029 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 2030 const int** offset, /**< pointer to store pointer that stores the offsets to each rows sparsity pattern in col, can be NULL */ 2031 const int** col /**< pointer to store pointer that stores the indices of variables that appear in each row, offset[nconss] gives length of col, can be NULL */ 2032 ) 2033 { 2034 SCIP_NLPIORACLECONS* cons; 2035 SCIP_EXPRITER* it; 2036 SCIP_EXPR* expr; 2037 SCIP_Bool* nzflag; 2038 int nnz; 2039 int maxnnz; 2040 int i; 2041 int j; 2042 2043 assert(oracle != NULL); 2044 2045 SCIPdebugMessage("%p get jacobian sparsity\n", (void*)oracle); 2046 2047 if( oracle->jacoffsets != NULL ) 2048 { 2049 assert(oracle->jaccols != NULL); 2050 if( offset != NULL ) 2051 *offset = oracle->jacoffsets; 2052 if( col != NULL ) 2053 *col = oracle->jaccols; 2054 return SCIP_OKAY; 2055 } 2056 2057 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 2058 2059 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &oracle->jacoffsets, oracle->nconss + 1) ); 2060 2061 maxnnz = MIN(oracle->nvars, 10) * oracle->nconss; /* initial guess */ 2062 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &oracle->jaccols, maxnnz) ); 2063 2064 if( maxnnz == 0 ) 2065 { 2066 /* no variables */ 2067 BMSclearMemoryArray(oracle->jacoffsets, oracle->nconss + 1); 2068 if( offset != NULL ) 2069 *offset = oracle->jacoffsets; 2070 if( col != NULL ) 2071 *col = oracle->jaccols; 2072 2073 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 2074 2075 return SCIP_OKAY; 2076 } 2077 nnz = 0; 2078 2079 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nzflag, oracle->nvars) ); 2080 2081 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 2082 SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) ); 2083 2084 for( i = 0; i < oracle->nconss; ++i ) 2085 { 2086 oracle->jacoffsets[i] = nnz; 2087 2088 cons = oracle->conss[i]; 2089 assert(cons != NULL); 2090 2091 if( cons->expr == NULL ) 2092 { 2093 /* for a linear constraint, we can just copy the linear indices from the constraint into the sparsity pattern */ 2094 if( cons->nlinidxs > 0 ) 2095 { 2096 SCIP_CALL( ensureIntArraySize(scip, &oracle->jaccols, &maxnnz, nnz + cons->nlinidxs) ); 2097 BMScopyMemoryArray(&oracle->jaccols[nnz], cons->linidxs, cons->nlinidxs); 2098 nnz += cons->nlinidxs; 2099 } 2100 continue; 2101 } 2102 2103 /* check which variables appear in constraint i 2104 * @todo this could be done faster for very sparse constraint by assembling all appearing variables, sorting, and removing duplicates 2105 */ 2106 BMSclearMemoryArray(nzflag, oracle->nvars); 2107 2108 for( j = 0; j < cons->nlinidxs; ++j ) 2109 nzflag[cons->linidxs[j]] = TRUE; 2110 2111 for( expr = SCIPexpriterRestartDFS(it, cons->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) 2112 if( SCIPisExprVaridx(scip, expr) ) 2113 { 2114 assert(SCIPgetIndexExprVaridx(expr) < oracle->nvars); 2115 nzflag[SCIPgetIndexExprVaridx(expr)] = TRUE; 2116 } 2117 2118 /* store variables indices in jaccols */ 2119 for( j = 0; j < oracle->nvars; ++j ) 2120 { 2121 if( nzflag[j] == FALSE ) 2122 continue; 2123 2124 SCIP_CALL( ensureIntArraySize(scip, &oracle->jaccols, &maxnnz, nnz + 1) ); 2125 oracle->jaccols[nnz] = j; 2126 ++nnz; 2127 } 2128 } 2129 2130 SCIPfreeExpriter(&it); 2131 2132 oracle->jacoffsets[oracle->nconss] = nnz; 2133 2134 /* shrink jaccols array to nnz */ 2135 if( nnz < maxnnz ) 2136 { 2137 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &oracle->jaccols, maxnnz, nnz) ); 2138 } 2139 2140 SCIPfreeBlockMemoryArray(scip, &nzflag, oracle->nvars); 2141 2142 if( offset != NULL ) 2143 *offset = oracle->jacoffsets; 2144 if( col != NULL ) 2145 *col = oracle->jaccols; 2146 2147 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 2148 2149 return SCIP_OKAY; 2150 } 2151 2152 /** evaluates the Jacobian matrix in a given point 2153 * 2154 * The values in the Jacobian matrix are returned in the same order as specified by the offset and col arrays obtained by SCIPnlpiOracleGetJacobianSparsity(). 2155 * The user need to call SCIPnlpiOracleGetJacobianSparsity() at least ones before using this function. 2156 * 2157 * @return SCIP_INVALIDDATA, if the Jacobian could not be evaluated (domain error, etc.) 2158 */ 2159 SCIP_RETCODE SCIPnlpiOracleEvalJacobian( 2160 SCIP* scip, /**< SCIP data structure */ 2161 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 2162 const SCIP_Real* x, /**< point where to evaluate */ 2163 SCIP_Bool isnewx, /**< has the point x changed since the last call to some evaluation function? */ 2164 SCIP_Real* convals, /**< pointer to store constraint values, can be NULL */ 2165 SCIP_Real* jacobi /**< pointer to store sparse jacobian values */ 2166 ) 2167 { 2168 SCIP_NLPIORACLECONS* cons; 2169 SCIP_RETCODE retcode; 2170 SCIP_Real* grad; 2171 SCIP_Real nlval; 2172 int i; 2173 int j; 2174 int k; 2175 int l; 2176 2177 SCIPdebugMessage("%p eval jacobian\n", (void*)oracle); 2178 2179 assert(oracle != NULL); 2180 assert(jacobi != NULL); 2181 2182 assert(oracle->jacoffsets != NULL); 2183 assert(oracle->jaccols != NULL); 2184 2185 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 2186 2187 SCIP_CALL( SCIPallocCleanBufferArray(scip, &grad, oracle->nvars) ); 2188 2189 retcode = SCIP_OKAY; 2190 2191 j = oracle->jacoffsets[0]; /* TODO isn't oracle->jacoffsets[0] == 0 and thus always j == k ? */ 2192 k = 0; 2193 for( i = 0; i < oracle->nconss; ++i ) 2194 { 2195 cons = oracle->conss[i]; 2196 assert(cons != NULL); 2197 2198 if( cons->expr == NULL ) 2199 { 2200 if( convals != NULL ) 2201 convals[i] = 0.0; 2202 2203 /* for a linear constraint, we can just copy the linear coefs from the constraint into the jacobian */ 2204 if( cons->nlinidxs > 0 ) 2205 { 2206 BMScopyMemoryArray(&jacobi[k], cons->lincoefs, cons->nlinidxs); 2207 j += cons->nlinidxs; 2208 k += cons->nlinidxs; 2209 if( convals != NULL ) 2210 for( l = 0; l < cons->nlinidxs; ++l ) 2211 convals[i] += cons->lincoefs[l] * x[cons->linidxs[l]]; 2212 } 2213 assert(j == oracle->jacoffsets[i+1]); 2214 continue; 2215 } 2216 2217 /* eval grad for nonlinear and add to jacobi */ 2218 SCIPdebugMsg(scip, "eval gradient of "); 2219 SCIPdebug( if( isnewx ) {printf("\nx ="); for( l = 0; l < oracle->nvars; ++l) printf(" %g", x[l]); printf("\n");} ) 2220 2221 SCIP_CALL( SCIPexprintGrad(scip, oracle->exprinterpreter, cons->expr, cons->exprintdata, (SCIP_Real*)x, isnewx, &nlval, grad) ); 2222 2223 SCIPdebug( printf("g ="); for( l = oracle->jacoffsets[i]; l < oracle->jacoffsets[i+1]; ++l) printf(" %g", grad[oracle->jaccols[l]]); printf("\n"); ) 2224 2225 if( !SCIPisFinite(nlval) || SCIPisInfinity(scip, ABS(nlval)) ) 2226 { 2227 SCIPdebugMessage("gradient evaluation yield invalid function value %g\n", nlval); 2228 retcode = SCIP_INVALIDDATA; /* indicate that the function could not be evaluated at given point */ 2229 goto TERMINATE; 2230 } 2231 if( convals != NULL ) 2232 convals[i] = nlval; 2233 2234 /* add linear part to grad */ 2235 for( l = 0; l < cons->nlinidxs; ++l ) 2236 { 2237 if( convals != NULL ) 2238 convals[i] += cons->lincoefs[l] * x[cons->linidxs[l]]; 2239 /* if grad[cons->linidxs[l]] is not finite, then adding a finite value doesn't change that, so don't check that here */ 2240 grad[cons->linidxs[l]] += cons->lincoefs[l]; 2241 } 2242 2243 /* store complete gradient (linear + nonlinear) in jacobi 2244 * use the already evaluated sparsity pattern to pick only elements from grad that could have been set 2245 */ 2246 assert(j == oracle->jacoffsets[i]); 2247 for( ; j < oracle->jacoffsets[i+1]; ++j ) 2248 { 2249 if( !SCIPisFinite(grad[oracle->jaccols[j]]) ) 2250 { 2251 SCIPdebugMessage("gradient evaluation yield invalid gradient value %g\n", grad[l]); 2252 retcode = SCIP_INVALIDDATA; /* indicate that the function could not be evaluated at given point */ 2253 goto TERMINATE; 2254 } 2255 jacobi[k++] = grad[oracle->jaccols[j]]; 2256 /* reset to 0 for next constraint */ 2257 grad[oracle->jaccols[j]] = 0.0; 2258 } 2259 2260 #ifndef NDEBUG 2261 /* check that exprint really wrote only into expected elements of grad 2262 * TODO remove after some testing for better performance of debug runs */ 2263 for( l = 0; l < oracle->nvars; ++l ) 2264 assert(grad[l] == 0.0); 2265 #endif 2266 } 2267 2268 TERMINATE: 2269 /* if there was an eval error, then we may have interrupted before cleaning up the grad buffer */ 2270 if( retcode == SCIP_INVALIDDATA ) 2271 BMSclearMemoryArray(grad, oracle->nvars); 2272 2273 SCIPfreeCleanBufferArray(scip, &grad); 2274 2275 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 2276 2277 return retcode; 2278 } 2279 2280 /** gets sparsity pattern of the Hessian matrix of the Lagrangian 2281 * 2282 * Note that internal data is returned in *offset and *col, thus the user must not to allocate memory there. 2283 * Adding or deleting variables, objective, or constraints may destroy the sparsity structure and make another call to this function necessary. 2284 * Only elements of the lower left triangle and the diagonal are counted. 2285 */ 2286 SCIP_RETCODE SCIPnlpiOracleGetHessianLagSparsity( 2287 SCIP* scip, /**< SCIP data structure */ 2288 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 2289 const int** offset, /**< pointer to store pointer that stores the offsets to each rows sparsity pattern in col, can be NULL */ 2290 const int** col /**< pointer to store pointer that stores the indices of variables that appear in each row, offset[nconss] gives length of col, can be NULL */ 2291 ) 2292 { 2293 int** colnz; /* nonzeros in Hessian corresponding to one column */ 2294 int* collen; /* collen[i] is length of array colnz[i] */ 2295 int* colnnz; /* colnnz[i] is number of entries in colnz[i] (<= collen[i]) */ 2296 int nnz; 2297 int i; 2298 int j; 2299 int cnt; 2300 2301 assert(oracle != NULL); 2302 2303 SCIPdebugMessage("%p get hessian lag sparsity\n", (void*)oracle); 2304 2305 if( oracle->heslagoffsets != NULL ) 2306 { 2307 assert(oracle->heslagcols != NULL); 2308 if( offset != NULL ) 2309 *offset = oracle->heslagoffsets; 2310 if( col != NULL ) 2311 *col = oracle->heslagcols; 2312 return SCIP_OKAY; 2313 } 2314 2315 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 2316 2317 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &oracle->heslagoffsets, oracle->nvars + 1) ); 2318 2319 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &colnz, oracle->nvars) ); 2320 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &collen, oracle->nvars) ); 2321 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &colnnz, oracle->nvars) ); 2322 BMSclearMemoryArray(colnz, oracle->nvars); 2323 BMSclearMemoryArray(collen, oracle->nvars); 2324 BMSclearMemoryArray(colnnz, oracle->nvars); 2325 nnz = 0; 2326 2327 if( oracle->objective->expr != NULL ) 2328 { 2329 SCIP_CALL( hessLagSparsitySetNzFlagForExpr(scip, oracle, colnz, collen, colnnz, &nnz, oracle->objective->expr, oracle->objective->exprintdata, oracle->nvars) ); 2330 } 2331 2332 for( i = 0; i < oracle->nconss; ++i ) 2333 { 2334 if( oracle->conss[i]->expr != NULL ) 2335 { 2336 SCIP_CALL( hessLagSparsitySetNzFlagForExpr(scip, oracle, colnz, collen, colnnz, &nnz, oracle->conss[i]->expr, oracle->conss[i]->exprintdata, oracle->nvars) ); 2337 } 2338 } 2339 2340 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &oracle->heslagcols, nnz) ); 2341 2342 /* set hessian sparsity from colnz, colnnz */ 2343 cnt = 0; 2344 for( i = 0; i < oracle->nvars; ++i ) 2345 { 2346 oracle->heslagoffsets[i] = cnt; 2347 for( j = 0; j < colnnz[i]; ++j ) 2348 { 2349 assert(cnt < nnz); 2350 oracle->heslagcols[cnt++] = colnz[i][j]; 2351 } 2352 SCIPfreeBlockMemoryArrayNull(scip, &colnz[i], collen[i]); 2353 collen[i] = 0; 2354 } 2355 oracle->heslagoffsets[oracle->nvars] = cnt; 2356 assert(cnt == nnz); 2357 2358 SCIPfreeBlockMemoryArray(scip, &colnz, oracle->nvars); 2359 SCIPfreeBlockMemoryArray(scip, &colnnz, oracle->nvars); 2360 SCIPfreeBlockMemoryArray(scip, &collen, oracle->nvars); 2361 2362 if( offset != NULL ) 2363 *offset = oracle->heslagoffsets; 2364 if( col != NULL ) 2365 *col = oracle->heslagcols; 2366 2367 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 2368 2369 return SCIP_OKAY; 2370 } 2371 2372 /** evaluates the Hessian matrix of the Lagrangian in a given point 2373 * 2374 * The values in the Hessian matrix are returned in the same order as specified by the offset and col arrays obtained by SCIPnlpiOracleGetHessianLagSparsity(). 2375 * The user must call SCIPnlpiOracleGetHessianLagSparsity() at least ones before using this function. 2376 * Only elements of the lower left triangle and the diagonal are computed. 2377 * 2378 * @return SCIP_INVALIDDATA, if the Hessian could not be evaluated (domain error, etc.) 2379 */ 2380 SCIP_RETCODE SCIPnlpiOracleEvalHessianLag( 2381 SCIP* scip, /**< SCIP data structure */ 2382 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 2383 const SCIP_Real* x, /**< point where to evaluate */ 2384 SCIP_Bool isnewx_obj, /**< has the point x changed since the last call to an objective evaluation function? */ 2385 SCIP_Bool isnewx_cons, /**< has the point x changed since the last call to the constraint evaluation function? */ 2386 SCIP_Real objfactor, /**< weight for objective function */ 2387 const SCIP_Real* lambda, /**< weights (Lagrangian multipliers) for the constraints */ 2388 SCIP_Real* hessian /**< pointer to store sparse hessian values */ 2389 ) 2390 { /*lint --e{715}*/ 2391 SCIP_RETCODE retcode = SCIP_OKAY; 2392 int i; 2393 2394 assert(oracle != NULL); 2395 assert(x != NULL); 2396 assert(lambda != NULL || oracle->nconss == 0); 2397 assert(hessian != NULL); 2398 2399 assert(oracle->heslagoffsets != NULL); 2400 assert(oracle->heslagcols != NULL); 2401 2402 SCIPdebugMessage("%p eval hessian lag\n", (void*)oracle); 2403 2404 SCIP_CALL( SCIPstartClock(scip, oracle->evalclock) ); 2405 2406 BMSclearMemoryArray(hessian, oracle->heslagoffsets[oracle->nvars]); 2407 2408 if( objfactor != 0.0 && oracle->objective->expr != NULL ) 2409 { 2410 retcode = hessLagAddExpr(scip, oracle, objfactor, x, isnewx_obj, oracle->objective->expr, oracle->objective->exprintdata, oracle->heslagoffsets, oracle->heslagcols, hessian); 2411 } 2412 2413 for( i = 0; i < oracle->nconss && retcode == SCIP_OKAY; ++i ) 2414 { 2415 assert( lambda != NULL ); /* for lint */ 2416 if( lambda[i] == 0.0 || oracle->conss[i]->expr == NULL ) 2417 continue; 2418 retcode = hessLagAddExpr(scip, oracle, lambda[i], x, isnewx_cons, oracle->conss[i]->expr, oracle->conss[i]->exprintdata, oracle->heslagoffsets, oracle->heslagcols, hessian); 2419 } 2420 2421 SCIP_CALL( SCIPstopClock(scip, oracle->evalclock) ); 2422 2423 return retcode; 2424 } 2425 2426 /** resets clock that measures evaluation time */ 2427 SCIP_RETCODE SCIPnlpiOracleResetEvalTime( 2428 SCIP* scip, /**< SCIP data structure */ 2429 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 2430 ) 2431 { 2432 assert(oracle != NULL); 2433 2434 SCIP_CALL( SCIPresetClock(scip, oracle->evalclock) ); 2435 2436 return SCIP_OKAY; 2437 } 2438 2439 /** gives time spend in evaluation since last reset of clock 2440 * 2441 * Gives 0 if the eval clock is disabled. 2442 */ 2443 SCIP_Real SCIPnlpiOracleGetEvalTime( 2444 SCIP* scip, /**< SCIP data structure */ 2445 SCIP_NLPIORACLE* oracle /**< pointer to NLPIORACLE data structure */ 2446 ) 2447 { 2448 assert(oracle != NULL); 2449 2450 return SCIPgetClockTime(scip, oracle->evalclock); 2451 } 2452 2453 /** prints the problem to a file. */ 2454 SCIP_RETCODE SCIPnlpiOraclePrintProblem( 2455 SCIP* scip, /**< SCIP data structure */ 2456 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 2457 FILE* file /**< file to print to, or NULL for standard output */ 2458 ) 2459 { /*lint --e{777} */ 2460 int i; 2461 SCIP_Real lhs; 2462 SCIP_Real rhs; 2463 2464 assert(oracle != NULL); 2465 2466 SCIPdebugMessage("%p print problem\n", (void*)oracle); 2467 2468 if( file == NULL ) 2469 file = stdout; 2470 2471 SCIPinfoMessage(scip, file, "NLPI Oracle %s: %d variables and %d constraints\n", oracle->name ? oracle->name : "", oracle->nvars, oracle->nconss); 2472 for( i = 0; i < oracle->nvars; ++i ) 2473 { 2474 if( oracle->varnames != NULL && oracle->varnames[i] != NULL ) 2475 SCIPinfoMessage(scip, file, "%10s (x%d)", oracle->varnames[i], i); /* give also name x%d as it will be by expression-print (printFunction) */ 2476 else 2477 SCIPinfoMessage(scip, file, "x%09d", i); 2478 SCIPinfoMessage(scip, file, ": [%8g, %8g]", oracle->varlbs[i], oracle->varubs[i]); 2479 SCIPinfoMessage(scip, file, "\t #linear: %d #nonlinear: %d\n", oracle->varlincount[i], oracle->varnlcount[i]); 2480 } 2481 2482 SCIPinfoMessage(scip, file, "objective: "); 2483 SCIP_CALL( printFunction(scip, oracle, file, oracle->objective, FALSE) ); 2484 if( oracle->objective->lhs != 0.0 ) 2485 SCIPinfoMessage(scip, file, "%+.15g", oracle->objective->lhs); 2486 SCIPinfoMessage(scip, file, "\n"); 2487 2488 for( i = 0; i < oracle->nconss; ++i ) 2489 { 2490 if( oracle->conss[i]->name != NULL ) 2491 SCIPinfoMessage(scip, file, "%10s", oracle->conss[i]->name); 2492 else 2493 SCIPinfoMessage(scip, file, "con%07d", i); 2494 2495 lhs = oracle->conss[i]->lhs; 2496 rhs = oracle->conss[i]->rhs; 2497 SCIPinfoMessage(scip, file, ": "); 2498 if( !SCIPisInfinity(scip, -lhs) && !SCIPisInfinity(scip, rhs) && lhs != rhs ) 2499 SCIPinfoMessage(scip, file, "%.15g <= ", lhs); 2500 2501 SCIP_CALL( printFunction(scip, oracle, file, oracle->conss[i], FALSE) ); 2502 2503 if( lhs == rhs ) 2504 SCIPinfoMessage(scip, file, " = %.15g", rhs); 2505 else if( !SCIPisInfinity(scip, rhs) ) 2506 SCIPinfoMessage(scip, file, " <= %.15g", rhs); 2507 else if( !SCIPisInfinity(scip, -lhs) ) 2508 SCIPinfoMessage(scip, file, " >= %.15g", lhs); 2509 2510 SCIPinfoMessage(scip, file, "\n"); 2511 } 2512 2513 return SCIP_OKAY; 2514 } 2515 2516 /** prints the problem to a file in GAMS format 2517 * 2518 * If there are variable (equation, resp.) names with more than 9 characters, then variable (equation, resp.) names are prefixed with an unique identifier. 2519 * This is to make it easier to identify variables solution output in the listing file. 2520 * Names with more than 64 characters are shorten to 64 letters due to GAMS limits. 2521 */ 2522 SCIP_RETCODE SCIPnlpiOraclePrintProblemGams( 2523 SCIP* scip, /**< SCIP data structure */ 2524 SCIP_NLPIORACLE* oracle, /**< pointer to NLPIORACLE data structure */ 2525 SCIP_Real* initval, /**< starting point values for variables or NULL */ 2526 FILE* file /**< file to print to, or NULL for standard output */ 2527 ) 2528 { /*lint --e{777} */ 2529 int i; 2530 int nllevel; /* level of nonlinearity of problem: linear = 0, quadratic, smooth nonlinear, nonsmooth */ 2531 static const char* nllevelname[4] = { "LP", "QCP", "NLP", "DNLP" }; 2532 char problemname[SCIP_MAXSTRLEN]; 2533 char namebuf[70]; 2534 SCIP_Bool havelongvarnames; 2535 SCIP_Bool havelongequnames; 2536 2537 SCIPdebugMessage("%p print problem gams\n", (void*)oracle); 2538 2539 assert(oracle != NULL); 2540 2541 if( file == NULL ) 2542 file = stdout; 2543 2544 nllevel = 0; 2545 2546 havelongvarnames = FALSE; 2547 for( i = 0; i < oracle->nvars; ++i ) 2548 if( oracle->varnames != NULL && oracle->varnames[i] != NULL && strlen(oracle->varnames[i]) > 9 ) 2549 { 2550 havelongvarnames = TRUE; 2551 break; 2552 } 2553 2554 havelongequnames = FALSE; 2555 for( i = 0; i < oracle->nconss; ++i ) 2556 if( oracle->conss[i]->name && strlen(oracle->conss[i]->name) > 9 ) 2557 { 2558 havelongequnames = TRUE; 2559 break; 2560 } 2561 2562 SCIPinfoMessage(scip, file, "$offlisting\n"); 2563 SCIPinfoMessage(scip, file, "$offdigit\n"); 2564 SCIPinfoMessage(scip, file, "* NLPI Oracle Problem %s\n", oracle->name ? oracle->name : ""); 2565 SCIPinfoMessage(scip, file, "Variables "); 2566 for( i = 0; i < oracle->nvars; ++i ) 2567 { 2568 printName(namebuf, oracle->varnames != NULL ? oracle->varnames[i] : NULL, i, 'x', NULL, havelongvarnames); 2569 SCIPinfoMessage(scip, file, "%s, ", namebuf); 2570 if( i % 10 == 9 ) 2571 SCIPinfoMessage(scip, file, "\n"); 2572 } 2573 SCIPinfoMessage(scip, file, "NLPIORACLEOBJVAR;\n\n"); 2574 for( i = 0; i < oracle->nvars; ++i ) 2575 { 2576 char* name; 2577 name = oracle->varnames != NULL ? oracle->varnames[i] : NULL; 2578 if( oracle->varlbs[i] == oracle->varubs[i] ) 2579 { 2580 printName(namebuf, name, i, 'x', NULL, havelongvarnames); 2581 SCIPinfoMessage(scip, file, "%s.fx = %.15g;\t", namebuf, oracle->varlbs[i]); 2582 } 2583 else 2584 { 2585 if( !SCIPisInfinity(scip, -oracle->varlbs[i]) ) 2586 { 2587 printName(namebuf, name, i, 'x', NULL, havelongvarnames); 2588 SCIPinfoMessage(scip, file, "%s.lo = %.15g;\t", namebuf, oracle->varlbs[i]); 2589 } 2590 if( !SCIPisInfinity(scip, oracle->varubs[i]) ) 2591 { 2592 printName(namebuf, name, i, 'x', NULL, havelongvarnames); 2593 SCIPinfoMessage(scip, file, "%s.up = %.15g;\t", namebuf, oracle->varubs[i]); 2594 } 2595 } 2596 if( initval != NULL ) 2597 { 2598 printName(namebuf, name, i, 'x', NULL, havelongvarnames); 2599 SCIPinfoMessage(scip, file, "%s.l = %.15g;\t", namebuf, initval[i]); 2600 } 2601 SCIPinfoMessage(scip, file, "\n"); 2602 } 2603 SCIPinfoMessage(scip, file, "\n"); 2604 2605 SCIPinfoMessage(scip, file, "Equations "); 2606 for( i = 0; i < oracle->nconss; ++i ) 2607 { 2608 printName(namebuf, oracle->conss[i]->name, i, 'e', NULL, havelongequnames); 2609 SCIPinfoMessage(scip, file, "%s, ", namebuf); 2610 2611 if( !SCIPisInfinity(scip, -oracle->conss[i]->lhs) && !SCIPisInfinity(scip, oracle->conss[i]->rhs) && oracle->conss[i]->lhs != oracle->conss[i]->rhs ) 2612 { 2613 /* ranged row: add second constraint */ 2614 printName(namebuf, oracle->conss[i]->name, i, 'e', "_RNG", havelongequnames); 2615 SCIPinfoMessage(scip, file, "%s, ", namebuf); 2616 } 2617 if( i % 10 == 9 ) 2618 SCIPinfoMessage(scip, file, "\n"); 2619 } 2620 SCIPinfoMessage(scip, file, "NLPIORACLEOBJ;\n\n"); 2621 2622 SCIPinfoMessage(scip, file, "NLPIORACLEOBJ.. NLPIORACLEOBJVAR =E= "); 2623 SCIP_CALL( printFunction(scip, oracle, file, oracle->objective, havelongvarnames) ); 2624 if( oracle->objective->lhs != 0.0 ) 2625 SCIPinfoMessage(scip, file, "%+.15g", oracle->objective->lhs); 2626 SCIPinfoMessage(scip, file, ";\n"); 2627 2628 for( i = 0; i < oracle->nconss; ++i ) 2629 { 2630 SCIP_Real lhs; 2631 SCIP_Real rhs; 2632 2633 printName(namebuf, oracle->conss[i]->name, i, 'e', NULL, havelongequnames); 2634 SCIPinfoMessage(scip, file, "%s.. ", namebuf); 2635 2636 SCIP_CALL( printFunction(scip, oracle, file, oracle->conss[i], havelongvarnames) ); 2637 2638 lhs = oracle->conss[i]->lhs; 2639 rhs = oracle->conss[i]->rhs; 2640 2641 if( lhs == rhs ) 2642 SCIPinfoMessage(scip, file, " =E= %.15g", rhs); 2643 else if( !SCIPisInfinity(scip, rhs) ) 2644 SCIPinfoMessage(scip, file, " =L= %.15g", rhs); 2645 else if( !SCIPisInfinity(scip, -lhs) ) 2646 SCIPinfoMessage(scip, file, " =G= %.15g", lhs); 2647 else 2648 SCIPinfoMessage(scip, file, " =N= 0"); 2649 SCIPinfoMessage(scip, file, ";\n"); 2650 2651 if( !SCIPisInfinity(scip, lhs) && !SCIPisInfinity(scip, rhs) && lhs != rhs ) 2652 { 2653 printName(namebuf, oracle->conss[i]->name, i, 'e', "_RNG", havelongequnames); 2654 SCIPinfoMessage(scip, file, "%s.. ", namebuf); 2655 2656 SCIP_CALL( printFunction(scip, oracle, file, oracle->conss[i], havelongvarnames) ); 2657 2658 SCIPinfoMessage(scip, file, " =G= %.15g;\n", lhs); 2659 } 2660 2661 if( nllevel <= 1 && oracle->conss[i]->expr != NULL ) 2662 nllevel = 2; 2663 if( nllevel <= 2 && oracle->conss[i]->expr != NULL ) 2664 { 2665 SCIP_Bool nonsmooth; 2666 SCIP_CALL( exprIsNonSmooth(scip, oracle->conss[i]->expr, &nonsmooth) ); 2667 if( nonsmooth ) 2668 nllevel = 3; 2669 } 2670 } 2671 2672 (void) SCIPsnprintf(problemname, SCIP_MAXSTRLEN, "%s", oracle->name ? oracle->name : "m"); 2673 2674 SCIPinfoMessage(scip, file, "Model %s / all /;\n", problemname); 2675 SCIPinfoMessage(scip, file, "option limrow = 0;\n"); 2676 SCIPinfoMessage(scip, file, "option limcol = 0;\n"); 2677 SCIPinfoMessage(scip, file, "Solve %s minimizing NLPIORACLEOBJVAR using %s;\n", problemname, nllevelname[nllevel]); 2678 2679 return SCIP_OKAY; 2680 } 2681 2682 /**@} */ 2683