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