1 /* Copyright (C) GAMS Development and others 2011
2 * All Rights Reserved.
3 * This code is published under the Eclipse Public License.
4 *
5 * Author: Stefan Vigerske
6 */
7
8 /**@file reader_gmo.c
9 * @ingroup FILEREADERS
10 * @brief GMO file reader
11 * @author Stefan Vigerske
12 */
13
14 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
15
16 #include <assert.h>
17 #include <math.h>
18 #include <string.h>
19
20 /* dos compiler does not know PI */
21 #ifndef M_PI
22 #define M_PI 3.141592653589793238462643
23 #endif
24
25 #include "reader_gmo.h"
26
27 #include "gmomcc.h"
28 #include "gevmcc.h"
29 #include "gdxcc.h"
30 #include "optcc.h"
31 #include "GamsCompatibility.h"
32 #include "GamsNLinstr.h"
33
34 #include "scip/cons_linear.h"
35 #include "scip/cons_bounddisjunction.h"
36 #include "scip/cons_quadratic.h"
37 #include "scip/cons_nonlinear.h"
38 #include "scip/cons_indicator.h"
39 #include "scip/cons_sos1.h"
40 #include "scip/cons_sos2.h"
41 #include "scip/heur_subnlp.h"
42 #include "scip/dialog_default.h"
43 #include "nlpi/struct_expr.h"
44
45 #define READER_NAME "gmoreader"
46 #define READER_DESC "Gams Control file reader (using GMO API)"
47 #define READER_EXTENSION "dat"
48
49
50 /*
51 * Data structures
52 */
53
54 /** data for gmo reader */
55 struct SCIP_ReaderData
56 {
57 gmoHandle_t gmo; /**< GAMS model object */
58 gevHandle_t gev; /**< GAMS environment */
59 int mipstart; /**< how to handle initial variable levels */
60 char* indicatorfile; /**< name of GAMS options file that contains definitions on indicators */
61 };
62
63 /** problem data */
64 struct SCIP_ProbData
65 {
66 int nvars; /**< number of variables in vars array */
67 SCIP_VAR** vars; /**< SCIP variables as corresponding to GMO variables */
68 SCIP_VAR* objvar; /**< SCIP variable used to model objective function */
69 SCIP_VAR* objconst; /**< SCIP variable used to model objective constant */
70 };
71
72 /*
73 * Callback methods of probdata
74 */
75
76 /** frees user data of original problem (called when the original problem is freed)
77 */
78 static
SCIP_DECL_PROBDELORIG(probdataDelOrigGmo)79 SCIP_DECL_PROBDELORIG(probdataDelOrigGmo)
80 {
81 int i;
82
83 assert((*probdata)->vars != NULL || (*probdata)->nvars > 0);
84
85 for( i = 0; i < (*probdata)->nvars; ++i )
86 {
87 SCIP_CALL( SCIPreleaseVar(scip, &(*probdata)->vars[i]) );
88 }
89 SCIPfreeMemoryArray(scip, &(*probdata)->vars);
90
91 if( (*probdata)->objvar != NULL )
92 {
93 SCIP_CALL( SCIPreleaseVar(scip, &(*probdata)->objvar) );
94 }
95
96 if( (*probdata)->objconst != NULL )
97 {
98 SCIP_CALL( SCIPreleaseVar(scip, &(*probdata)->objconst) );
99 }
100
101 SCIPfreeMemory(scip, probdata);
102
103 return SCIP_OKAY;
104 }
105
106 /** creates user data of transformed problem by transforming the original user problem data
107 * (called after problem was transformed)
108 */
109 #if 0
110 static
111 SCIP_DECL_PROBTRANS(probdataTransGmo)
112 {
113
114 }
115 #else
116 #define probdataTransGmo NULL
117 #endif
118
119 /** frees user data of transformed problem (called when the transformed problem is freed)
120 */
121 #if 0
122 static
123 SCIP_DECL_PROBDELTRANS(probdataDelTransGmo)
124 {
125 return SCIP_OKAY;
126 }
127 #else
128 #define probdataDelTransGmo NULL
129 #endif
130
131
132 /** solving process initialization method of transformed data (called before the branch and bound process begins)
133 */
134 #if 0
135 static
136 SCIP_DECL_PROBINITSOL(probdataInitSolGmo)
137 {
138
139 }
140 #else
141 #define probdataInitSolGmo NULL
142 #endif
143
144 /** solving process deinitialization method of transformed data (called before the branch and bound data is freed)
145 */
146 #if 0
147 static
148 SCIP_DECL_PROBEXITSOL(probdataExitSolGmo)
149 {
150 return SCIP_OKAY;
151 }
152 #else
153 #define probdataExitSolGmo NULL
154 #endif
155
156 /** copies user data of source SCIP for the target SCIP
157 */
158 #if 0
159 static
160 SCIP_DECL_PROBCOPY(probdataCopyGmo)
161 {
162
163 }
164 #else
165 #define probdataCopyGmo NULL
166 #endif
167
168 /*
169 * Local methods of reader
170 */
171
172
173 /** ensures that an array of variables has at least a given length */
174 static
ensureVarsSize(SCIP * scip,SCIP_VAR *** vars,int * varssize,int nvars)175 SCIP_RETCODE ensureVarsSize(
176 SCIP* scip, /**< SCIP data structure */
177 SCIP_VAR*** vars, /**< pointer to variables array */
178 int* varssize, /**< pointer to length of variables array */
179 int nvars /**< desired minimal length of array */
180 )
181 {
182 assert(scip != NULL);
183 assert(vars != NULL);
184 assert(varssize != NULL);
185
186 if( nvars < *varssize )
187 return SCIP_OKAY;
188
189 *varssize = SCIPcalcMemGrowSize(scip, nvars);
190 SCIP_CALL( SCIPreallocBufferArray(scip, vars, *varssize) );
191 assert(nvars <= *varssize);
192
193 return SCIP_OKAY;
194 }
195
196 /** adds new children to a linear expression */
197 static
exprLinearAdd(BMS_BLKMEM * blkmem,SCIP_EXPR * expr,int nchildren,SCIP_Real * coefs,SCIP_EXPR ** children,SCIP_Real constant)198 SCIP_RETCODE exprLinearAdd(
199 BMS_BLKMEM* blkmem, /**< block memory */
200 SCIP_EXPR* expr, /**< pointer to store resulting expression */
201 int nchildren, /**< number of children to add */
202 SCIP_Real* coefs, /**< coefficients of additional children */
203 SCIP_EXPR** children, /**< additional children expressions */
204 SCIP_Real constant /**< constant to add */
205 )
206 {
207 SCIP_Real* data;
208
209 assert(blkmem != NULL);
210 assert(expr != NULL);
211 assert(SCIPexprGetOperator(expr) == SCIP_EXPR_LINEAR);
212 assert(nchildren >= 0);
213 assert(coefs != NULL || nchildren == 0);
214 assert(children != NULL || nchildren == 0);
215
216 data = (SCIP_Real*)SCIPexprGetOpData(expr);
217 assert(data != NULL);
218
219 /* handle simple case of adding a constant */
220 if( nchildren == 0 )
221 {
222 data[SCIPexprGetNChildren(expr)] += constant;
223
224 return SCIP_OKAY;
225 }
226
227 /* add new children to expr's children array */
228 SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, expr->nchildren + nchildren) );
229 BMScopyMemoryArray(&expr->children[expr->nchildren], children, nchildren); /*lint !e866*/
230
231 /* add constant and new coefs to expr's data array */
232 SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &data, expr->nchildren + 1, expr->nchildren + nchildren + 1) );
233 data[expr->nchildren + nchildren] = data[expr->nchildren] + constant;
234 BMScopyMemoryArray(&data[expr->nchildren], coefs, nchildren); /*lint !e866*/
235 expr->data.data = (void*)data;
236
237 expr->nchildren += nchildren;
238
239 return SCIP_OKAY;
240 }
241
242 /** frees a linear expression, but not its children */
243 static
exprLinearFree(BMS_BLKMEM * blkmem,SCIP_EXPR ** expr)244 void exprLinearFree(
245 BMS_BLKMEM* blkmem, /**< block memory */
246 SCIP_EXPR** expr /**< linear expression to free */
247 )
248 {
249 assert(blkmem != NULL);
250 assert(expr != NULL);
251 assert(SCIPexprGetOperator(*expr) == SCIP_EXPR_LINEAR);
252
253 BMSfreeBlockMemoryArray(blkmem, (SCIP_Real**)&(*expr)->data.data, (*expr)->nchildren + 1); /*lint !e866*/
254 BMSfreeBlockMemoryArray(blkmem, &(*expr)->children, (*expr)->nchildren);
255
256 BMSfreeBlockMemory(blkmem, expr);
257 }
258
259 /** creates an expression from the addition of two given expression, with coefficients, and a constant */
260 static
exprAdd(BMS_BLKMEM * blkmem,SCIP_EXPR ** expr,SCIP_Real coef1,SCIP_EXPR * term1,SCIP_Real coef2,SCIP_EXPR * term2,SCIP_Real constant)261 SCIP_RETCODE exprAdd(
262 BMS_BLKMEM* blkmem, /**< block memory */
263 SCIP_EXPR** expr, /**< pointer to store resulting expression */
264 SCIP_Real coef1, /**< coefficient of first term */
265 SCIP_EXPR* term1, /**< expression of first term */
266 SCIP_Real coef2, /**< coefficient of second term */
267 SCIP_EXPR* term2, /**< expression of second term */
268 SCIP_Real constant /**< constant term to add */
269 )
270 {
271 assert(blkmem != NULL);
272 assert(expr != NULL);
273
274 if( term1 != NULL && SCIPexprGetOperator(term1) == SCIP_EXPR_CONST )
275 {
276 constant += coef1 * SCIPexprGetOpReal(term1);
277 SCIPexprFreeDeep(blkmem, &term1);
278 }
279
280 if( term2 != NULL && SCIPexprGetOperator(term2) == SCIP_EXPR_CONST )
281 {
282 constant += coef2 * SCIPexprGetOpReal(term2);
283 SCIPexprFreeDeep(blkmem, &term2);
284 }
285
286 if( term1 == NULL && term2 == NULL )
287 {
288 SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_CONST, constant) );
289 return SCIP_OKAY;
290 }
291
292 if( term1 != NULL && SCIPexprGetOperator(term1) == SCIP_EXPR_LINEAR && coef1 != 1.0 )
293 {
294 /* multiply coefficients and constant of linear expression term1 by coef1 */
295 SCIP_Real* coefs;
296 int i;
297
298 coefs = SCIPexprGetLinearCoefs(term1);
299 assert(coefs != NULL);
300
301 for( i = 0; i < SCIPexprGetNChildren(term1); ++i )
302 coefs[i] *= coef1;
303
304 SCIP_CALL( exprLinearAdd(blkmem, term1, 0, NULL, NULL, (coef1-1.0) * SCIPexprGetLinearConstant(term1)) );
305
306 coef1 = 1.0;
307 }
308
309 if( term2 != NULL && SCIPexprGetOperator(term2) == SCIP_EXPR_LINEAR && coef2 != 1.0 )
310 {
311 /* multiply coefficients and constant of linear expression term2 by coef2 */
312 SCIP_Real* coefs;
313 int i;
314
315 coefs = SCIPexprGetLinearCoefs(term2);
316 assert(coefs != NULL);
317
318 for( i = 0; i < SCIPexprGetNChildren(term2); ++i )
319 coefs[i] *= coef2;
320
321 SCIP_CALL( exprLinearAdd(blkmem, term2, 0, NULL, NULL, (coef2-1.0) * SCIPexprGetLinearConstant(term2)) );
322
323 coef2 = 1.0;
324 }
325
326 if( term1 == NULL || term2 == NULL )
327 {
328 if( term1 == NULL )
329 {
330 term1 = term2;
331 coef1 = coef2;
332 /* term2 = NULL; */
333 }
334 if( constant != 0.0 || coef1 != 1.0 )
335 {
336 if( SCIPexprGetOperator(term1) == SCIP_EXPR_LINEAR )
337 {
338 assert(coef1 == 1.0);
339
340 /* add constant to existing linear expression */
341 SCIP_CALL( exprLinearAdd(blkmem, term1, 0, NULL, NULL, constant) );
342 *expr = term1;
343 }
344 else
345 {
346 /* create new linear expression for coef1 * term1 + constant */
347 SCIP_CALL( SCIPexprCreateLinear(blkmem, expr, 1, &term1, &coef1, constant) );
348 }
349 }
350 else
351 {
352 assert(constant == 0.0);
353 assert(coef1 == 1.0);
354 *expr = term1;
355 }
356
357 return SCIP_OKAY;
358 }
359
360 if( SCIPexprGetOperator(term1) == SCIP_EXPR_LINEAR && SCIPexprGetOperator(term2) == SCIP_EXPR_LINEAR )
361 {
362 assert(coef1 == 1.0);
363 assert(coef2 == 1.0);
364
365 SCIP_CALL( exprLinearAdd(blkmem, term1, SCIPexprGetNChildren(term2), SCIPexprGetLinearCoefs(term2), SCIPexprGetChildren(term2), SCIPexprGetLinearConstant(term2) + constant) );
366 exprLinearFree(blkmem, &term2);
367
368 *expr = term1;
369
370 return SCIP_OKAY;
371 }
372
373 if( SCIPexprGetOperator(term2) == SCIP_EXPR_LINEAR )
374 {
375 /* if only term2 is linear, then swap */
376 SCIP_EXPR* tmp;
377
378 tmp = term2;
379 assert(coef2 == 1.0);
380
381 term2 = term1;
382 coef2 = coef1;
383 term1 = tmp;
384 coef1 = 1.0;
385 }
386
387 if( SCIPexprGetOperator(term1) == SCIP_EXPR_LINEAR )
388 {
389 /* add coef2*term2 as extra child to linear expression term1 */
390 assert(coef1 == 1.0);
391
392 SCIP_CALL( exprLinearAdd(blkmem, term1, 1, &coef2, &term2, constant) );
393 *expr = term1;
394
395 return SCIP_OKAY;
396 }
397
398 /* both terms are not linear, then create new linear term for sum */
399 {
400 SCIP_Real coefs[2];
401 SCIP_EXPR* children[2];
402
403 coefs[0] = coef1;
404 coefs[1] = coef2;
405 children[0] = term1;
406 children[1] = term2;
407
408 SCIP_CALL( SCIPexprCreateLinear(blkmem, expr, 2, children, coefs, constant) );
409 }
410
411 return SCIP_OKAY;
412 }
413
414 /** creates an expression tree from given GAMS nonlinear instructions */
415 static
makeExprtree(SCIP * scip,gmoHandle_t gmo,int codelen,int * opcodes,int * fields,SCIP_Real * constants,SCIP_EXPRTREE ** exprtree)416 SCIP_RETCODE makeExprtree(
417 SCIP* scip, /**< SCIP data structure */
418 gmoHandle_t gmo, /**< GAMS Model Object */
419 int codelen,
420 int* opcodes,
421 int* fields,
422 SCIP_Real* constants,
423 SCIP_EXPRTREE** exprtree /**< buffer where to store expression tree */
424 )
425 {
426 SCIP_PROBDATA* probdata;
427 BMS_BLKMEM* blkmem;
428 SCIP_HASHMAP* var2idx;
429 SCIP_EXPR** stack;
430 int stackpos;
431 int stacksize;
432 SCIP_VAR** vars;
433 int nvars;
434 int varssize;
435 int pos;
436 SCIP_EXPR* e;
437 SCIP_EXPR* term1;
438 SCIP_EXPR* term2;
439 GamsOpCode opcode;
440 int address;
441 int varidx;
442 int nargs;
443 int rc;
444
445 assert(scip != NULL);
446 assert(opcodes != NULL);
447 assert(fields != NULL);
448 assert(constants != NULL);
449 assert(exprtree != NULL);
450
451 probdata = SCIPgetProbData(scip);
452 assert(probdata != NULL);
453 assert(probdata->vars != NULL);
454
455 blkmem = SCIPblkmem(scip);
456
457 stackpos = 0;
458 stacksize = 20;
459 SCIP_CALL( SCIPallocBufferArray(scip, &stack, stacksize) );
460
461 nvars = 0;
462 varssize = 10;
463 SCIP_CALL( SCIPallocBufferArray(scip, &vars, varssize) );
464 SCIP_CALL( SCIPhashmapCreate(&var2idx, blkmem, SCIPgetNVars(scip)) );
465
466 nargs = -1;
467 rc = SCIP_OKAY;
468
469 for( pos = 0; pos < codelen; ++pos )
470 {
471 opcode = (GamsOpCode)opcodes[pos];
472 address = fields[pos]-1;
473
474 SCIPdebugMessage("%s: ", GamsOpCodeName[opcode]);
475
476 e = NULL;
477
478 switch( opcode )
479 {
480 case nlNoOp: /* no operation */
481 case nlStore: /* store row */
482 case nlHeader:
483 {
484 SCIPdebugPrintf("ignored\n");
485 break;
486 }
487
488 case nlPushV: /* push variable */
489 {
490 address = gmoGetjSolver(gmo, address);
491 SCIPdebugPrintf("push variable %d = <%s>\n", address, SCIPvarGetName(probdata->vars[address]));
492
493 if( !SCIPhashmapExists(var2idx, probdata->vars[address]) )
494 {
495 /* add variable to list of variables */
496 SCIP_CALL( ensureVarsSize(scip, &vars, &varssize, nvars+1) );
497 assert(nvars < varssize);
498 vars[nvars] = probdata->vars[address];
499 varidx = nvars;
500 ++nvars;
501 SCIP_CALL( SCIPhashmapInsert(var2idx, (void*)vars[varidx], (void*)(size_t)varidx) );
502 }
503 else
504 {
505 varidx = (int)(size_t)SCIPhashmapGetImage(var2idx, (void*)probdata->vars[address]);
506 assert(varidx >= 0);
507 assert(varidx < nvars);
508 assert(vars[varidx] == probdata->vars[address]);
509 }
510 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_VARIDX, varidx) );
511 break;
512 }
513
514 case nlPushI: /* push constant */
515 {
516 SCIPdebugPrintf("push constant %g\n", constants[address]);
517 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_CONST, constants[address]) );
518 break;
519 }
520
521 case nlPushZero: /* push zero */
522 {
523 SCIPdebugPrintf("push constant zero\n");
524
525 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_CONST, 0.0) );
526 break;
527 }
528
529 case nlAdd : /* add */
530 {
531 SCIPdebugPrintf("add\n");
532 assert(stackpos >= 2);
533 term1 = stack[stackpos-1];
534 --stackpos;
535 term2 = stack[stackpos-1];
536 --stackpos;
537
538 SCIP_CALL( exprAdd(blkmem, &e, 1.0, term1, 1.0, term2, 0.0) );
539
540 break;
541 }
542
543 case nlAddV: /* add variable */
544 {
545 address = gmoGetjSolver(gmo, address);
546 SCIPdebugPrintf("add variable %d = <%s>\n", address, SCIPvarGetName(probdata->vars[address]));
547
548 assert(stackpos >= 1);
549 term1 = stack[stackpos-1];
550 --stackpos;
551
552 if( !SCIPhashmapExists(var2idx, probdata->vars[address]) )
553 {
554 /* add variable to list of variables */
555 SCIP_CALL( ensureVarsSize(scip, &vars, &varssize, nvars+1) );
556 assert(nvars < varssize);
557 vars[nvars] = probdata->vars[address];
558 varidx = nvars;
559 ++nvars;
560 SCIP_CALL( SCIPhashmapInsert(var2idx, (void*)vars[varidx], (void*)(size_t)varidx) );
561 }
562 else
563 {
564 varidx = (int)(size_t)SCIPhashmapGetImage(var2idx, (void*)probdata->vars[address]);
565 assert(varidx >= 0);
566 assert(varidx < nvars);
567 assert(vars[varidx] == probdata->vars[address]);
568 }
569
570 SCIP_CALL( SCIPexprCreate(blkmem, &term2, SCIP_EXPR_VARIDX, varidx) );
571 SCIP_CALL( exprAdd(blkmem, &e, 1.0, term1, 1.0, term2, 0.0) );
572
573 break;
574 }
575
576 case nlAddI: /* add immediate */
577 {
578 SCIPdebugPrintf("add constant %g\n", constants[address]);
579
580 assert(stackpos >= 1);
581 term1 = stack[stackpos-1];
582 --stackpos;
583
584 SCIP_CALL( exprAdd(blkmem, &e, 1.0, term1, 1.0, NULL, constants[address]) );
585
586 break;
587 }
588
589 case nlSub: /* substract */
590 {
591 SCIPdebugPrintf("substract\n");
592
593 assert(stackpos >= 2);
594 term1 = stack[stackpos-1];
595 --stackpos;
596 term2 = stack[stackpos-1];
597 --stackpos;
598
599 SCIP_CALL( exprAdd(blkmem, &e, 1.0, term2, -1.0, term1, 0.0) );
600
601 break;
602 }
603
604 case nlSubV: /* subtract variable */
605 {
606 address = gmoGetjSolver(gmo, address);
607 SCIPdebugPrintf("subtract variable %d = <%s>\n", address, SCIPvarGetName(probdata->vars[address]));
608
609 assert(stackpos >= 1);
610 term1 = stack[stackpos-1];
611 --stackpos;
612
613 if( !SCIPhashmapExists(var2idx, probdata->vars[address]) )
614 {
615 /* add variable to list of variables */
616 SCIP_CALL( ensureVarsSize(scip, &vars, &varssize, nvars+1) );
617 assert(nvars < varssize);
618 vars[nvars] = probdata->vars[address];
619 varidx = nvars;
620 ++nvars;
621 SCIP_CALL( SCIPhashmapInsert(var2idx, (void*)vars[varidx], (void*)(size_t)varidx) );
622 }
623 else
624 {
625 varidx = (int)(size_t)SCIPhashmapGetImage(var2idx, (void*)probdata->vars[address]);
626 assert(varidx >= 0);
627 assert(varidx < nvars);
628 assert(vars[varidx] == probdata->vars[address]);
629 }
630
631 SCIP_CALL( SCIPexprCreate(blkmem, &term2, SCIP_EXPR_VARIDX, varidx) );
632 SCIP_CALL( exprAdd(blkmem, &e, 1.0, term1, -1.0, term2, 0.0) );
633
634 break;
635 }
636
637 case nlSubI: /* subtract immediate */
638 {
639 SCIPdebugPrintf("subtract constant %g\n", constants[address]);
640
641 assert(stackpos >= 1);
642 term1 = stack[stackpos-1];
643 --stackpos;
644
645 SCIP_CALL( exprAdd(blkmem, &e, 1.0, term1, 1.0, NULL, -constants[address]) );
646
647 break;
648 }
649
650 case nlMul: /* multiply */
651 {
652 SCIPdebugPrintf("multiply\n");
653
654 assert(stackpos >= 2);
655 term1 = stack[stackpos-1];
656 --stackpos;
657 term2 = stack[stackpos-1];
658 --stackpos;
659
660 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MUL, term2, term1) );
661 break;
662 }
663
664 case nlMulV: /* multiply variable */
665 {
666 address = gmoGetjSolver(gmo, address);
667 SCIPdebugPrintf("multiply variable %d = <%s>\n", address, SCIPvarGetName(probdata->vars[address]));
668
669 assert(stackpos >= 1);
670 term1 = stack[stackpos-1];
671 --stackpos;
672
673 if( !SCIPhashmapExists(var2idx, probdata->vars[address]) )
674 {
675 /* add variable to list of variables */
676 SCIP_CALL( ensureVarsSize(scip, &vars, &varssize, nvars+1) );
677 assert(nvars < varssize);
678 vars[nvars] = probdata->vars[address];
679 varidx = nvars;
680 ++nvars;
681 SCIP_CALL( SCIPhashmapInsert(var2idx, (void*)vars[varidx], (void*)(size_t)varidx) );
682 }
683 else
684 {
685 varidx = (int)(size_t)SCIPhashmapGetImage(var2idx, (void*)probdata->vars[address]);
686 assert(varidx >= 0);
687 assert(varidx < nvars);
688 assert(vars[varidx] == probdata->vars[address]);
689 }
690
691 SCIP_CALL( SCIPexprCreate(blkmem, &term2, SCIP_EXPR_VARIDX, varidx) );
692 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MUL, term1, term2) );
693 break;
694 }
695
696 case nlMulI: /* multiply immediate */
697 {
698 SCIPdebugPrintf("multiply constant %g\n", constants[address]);
699
700 assert(stackpos >= 1);
701 term1 = stack[stackpos-1];
702 --stackpos;
703
704 SCIP_CALL( exprAdd(blkmem, &e, constants[address], term1, 1.0, NULL, 0.0) );
705
706 break;
707 }
708
709 case nlMulIAdd:
710 {
711 SCIPdebugPrintf("multiply constant %g and add\n", constants[address]);
712
713 assert(stackpos >= 2);
714 term1 = stack[stackpos-1];
715 --stackpos;
716 term2 = stack[stackpos-1];
717 --stackpos;
718
719 SCIP_CALL( exprAdd(blkmem, &e, constants[address], term1, 1.0, term2, 0.0) );
720
721 break;
722 }
723
724 #if 1
725 case nlDiv: /* divide */
726 {
727 SCIPdebugPrintf("divide\n");
728
729 assert(stackpos >= 2);
730 term1 = stack[stackpos-1];
731 --stackpos;
732 term2 = stack[stackpos-1];
733 --stackpos;
734
735 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_DIV, term2, term1) );
736 break;
737 }
738 #elif 1
739 case nlDiv: /* divide */
740 {
741 SCIP_EXPRDATA_MONOMIAL* monomial;
742 SCIP_EXPR* children[2];
743 double exps[2] = { -1.0, 1.0 };
744
745 SCIPdebugPrintf("divide\n");
746
747 assert(stackpos >= 2);
748 children[0] = stack[stackpos-1];
749 --stackpos;
750 children[1] = stack[stackpos-1];
751 --stackpos;
752
753 SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 2, NULL, exps) );
754 SCIP_CALL( SCIPexprCreatePolynomial(blkmem, &e, 2, children, 1, &monomial, 0.0, FALSE) );
755 break;
756 }
757 #else
758 case nlDiv: /* divide */
759 {
760 SCIPdebugPrintf("divide\n");
761
762 assert(stackpos >= 2);
763 term1 = stack[stackpos-1];
764 --stackpos;
765 term2 = stack[stackpos-1];
766 --stackpos;
767
768 SCIP_CALL( SCIPexprCreate(blkmem, &term1, SCIP_EXPR_INTPOWER, term1, -1) );
769
770 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MUL, term2, term1) );
771 break;
772 }
773 #endif
774
775 case nlDivV: /* divide variable */
776 {
777 address = gmoGetjSolver(gmo, address);
778 SCIPdebugPrintf("divide variable %d = <%s>\n", address, SCIPvarGetName(probdata->vars[address]));
779
780 assert(stackpos >= 1);
781 term1 = stack[stackpos-1];
782 --stackpos;
783
784 if( !SCIPhashmapExists(var2idx, probdata->vars[address]) )
785 {
786 /* add variable to list of variables */
787 SCIP_CALL( ensureVarsSize(scip, &vars, &varssize, nvars+1) );
788 assert(nvars < varssize);
789 vars[nvars] = probdata->vars[address];
790 varidx = nvars;
791 ++nvars;
792 SCIP_CALL( SCIPhashmapInsert(var2idx, (void*)vars[varidx], (void*)(size_t)varidx) );
793 }
794 else
795 {
796 varidx = (int)(size_t)SCIPhashmapGetImage(var2idx, (void*)probdata->vars[address]);
797 assert(varidx >= 0);
798 assert(varidx < nvars);
799 assert(vars[varidx] == probdata->vars[address]);
800 }
801
802 SCIP_CALL( SCIPexprCreate(blkmem, &term2, SCIP_EXPR_VARIDX, varidx) );
803 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_DIV, term1, term2) );
804 break;
805 }
806
807 case nlDivI: /* divide immediate */
808 {
809 SCIPdebugPrintf("divide constant %g\n", constants[address]);
810 assert(constants[address] != 0.0);
811
812 assert(stackpos >= 1);
813 term1 = stack[stackpos-1];
814 --stackpos;
815
816 SCIP_CALL( exprAdd(blkmem, &e, 1.0/constants[address], term1, 1.0, NULL, 0.0) );
817
818 break;
819 }
820
821 case nlUMin: /* unary minus */
822 {
823 SCIPdebugPrintf("negate\n");
824
825 assert(stackpos >= 1);
826 term1 = stack[stackpos-1];
827 --stackpos;
828
829 SCIP_CALL( exprAdd(blkmem, &e, -1.0, term1, 1.0, NULL, 0.0) );
830
831 break;
832 }
833
834 case nlUMinV: /* unary minus variable */
835 {
836 SCIP_Real minusone;
837
838 address = gmoGetjSolver(gmo, address);
839 SCIPdebugPrintf("push negated variable %d = <%s>\n", address, SCIPvarGetName(probdata->vars[address]));
840
841 if( !SCIPhashmapExists(var2idx, probdata->vars[address]) )
842 {
843 /* add variable to list of variables */
844 SCIP_CALL( ensureVarsSize(scip, &vars, &varssize, nvars+1) );
845 assert(nvars < varssize);
846 vars[nvars] = probdata->vars[address];
847 varidx = nvars;
848 ++nvars;
849 SCIP_CALL( SCIPhashmapInsert(var2idx, (void*)vars[varidx], (void*)(size_t)varidx) );
850 }
851 else
852 {
853 varidx = (int)(size_t)SCIPhashmapGetImage(var2idx, (void*)probdata->vars[address]);
854 assert(varidx >= 0);
855 assert(varidx < nvars);
856 assert(vars[varidx] == probdata->vars[address]);
857 }
858
859 minusone = -1.0;
860 SCIP_CALL( SCIPexprCreate(blkmem, &term1, SCIP_EXPR_VARIDX, varidx) );
861 SCIP_CALL( SCIPexprCreateLinear(blkmem, &e, 1, &term1, &minusone, 0.0) );
862
863 break;
864 }
865
866 case nlFuncArgN:
867 {
868 SCIPdebugPrintf("set number of arguments = %d\n", address);
869 nargs = address + 1; /* undo shift by one */
870 break;
871 }
872
873 case nlCallArg1:
874 nargs = 1;
875 /*lint -fallthrough*/
876
877 case nlCallArg2:
878 if( opcode == nlCallArg2 )
879 nargs = 2;
880 /*lint -fallthrough*/
881
882 case nlCallArgN:
883 {
884 GamsFuncCode func;
885
886 SCIPdebugPrintf("call function ");
887
888 func = (GamsFuncCode)(address+1); /* here the shift by one was not a good idea */
889
890 switch( func )
891 {
892 case fnmin:
893 {
894 SCIPdebugPrintf("min\n");
895
896 assert(stackpos >= 2);
897 term1 = stack[stackpos-1];
898 --stackpos;
899 term2 = stack[stackpos-1];
900 --stackpos;
901
902 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MIN, term1, term2) );
903 break;
904 }
905
906 case fnmax:
907 {
908 SCIPdebugPrintf("max\n");
909
910 assert(stackpos >= 2);
911 term1 = stack[stackpos-1];
912 --stackpos;
913 term2 = stack[stackpos-1];
914 --stackpos;
915
916 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MAX, term1, term2) );
917 break;
918 }
919
920 case fnsqr:
921 {
922 SCIPdebugPrintf("square\n");
923
924 assert(stackpos >= 1);
925 term1 = stack[stackpos-1];
926 --stackpos;
927
928 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_SQUARE, term1) );
929 break;
930 }
931
932 case fnexp:
933 {
934 SCIPdebugPrintf("exp\n");
935
936 assert(stackpos >= 1);
937 term1 = stack[stackpos-1];
938 --stackpos;
939
940 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_EXP, term1) );
941 break;
942 }
943
944 case fnlog:
945 {
946 SCIPdebugPrintf("log\n");
947
948 assert(stackpos >= 1);
949 term1 = stack[stackpos-1];
950 --stackpos;
951
952 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_LOG, term1) );
953 break;
954 }
955
956 case fnlog10:
957 {
958 SCIPdebugPrintf("log10 = ln * 1/ln(10)\n");
959
960 assert(stackpos >= 1);
961 term1 = stack[stackpos-1];
962 --stackpos;
963
964 SCIP_CALL( SCIPexprCreate(blkmem, &term2, SCIP_EXPR_LOG, term1) );
965 SCIP_CALL( SCIPexprCreate(blkmem, &term1, SCIP_EXPR_CONST, 1.0/log(10.0)) );
966 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MUL, term2, term1) );
967 break;
968 }
969
970 case fnlog2:
971 {
972 SCIPdebugPrintf("log2 = ln * 1/ln(2)\n");
973
974 assert(stackpos >= 1);
975 term1 = stack[stackpos-1];
976 --stackpos;
977
978 SCIP_CALL( SCIPexprCreate(blkmem, &term2, SCIP_EXPR_LOG, term1) );
979 SCIP_CALL( SCIPexprCreate(blkmem, &term1, SCIP_EXPR_CONST, 1.0/log(2.0)) );
980 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MUL, term2, term1) );
981 break;
982 }
983
984 case fnsqrt:
985 {
986 SCIPdebugPrintf("sqrt\n");
987
988 assert(stackpos >= 1);
989 term1 = stack[stackpos-1];
990 --stackpos;
991
992 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_SQRT, term1) );
993 break;
994 }
995 #if 0
996 case fncos:
997 {
998 SCIPdebugPrintf("cos\n");
999
1000 assert(stackpos >= 1);
1001 term1 = stack[stackpos-1];
1002 --stackpos;
1003
1004 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_COS, term1) );
1005 break;
1006 }
1007
1008 case fnsin:
1009 {
1010 SCIPdebugPrintf("sin\n");
1011
1012 assert(stackpos >= 1);
1013 term1 = stack[stackpos-1];
1014 --stackpos;
1015
1016 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_SIN, term1) );
1017 break;
1018 }
1019
1020 case fntan:
1021 {
1022 SCIPdebugPrintf("tan\n");
1023
1024 assert(stackpos >= 1);
1025 term1 = stack[stackpos-1];
1026 --stackpos;
1027
1028 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_TAN, term1) );
1029 break;
1030 }
1031 #endif
1032 case fnpower: /* x ^ y */
1033 case fnrpower: /* x ^ y */
1034 case fncvpower: /* constant ^ x */
1035 case fnvcpower: /* x ^ constant */
1036 {
1037 SCIPdebugPrintf("power\n");
1038
1039 assert(stackpos >= 2);
1040 term1 = stack[stackpos-1];
1041 --stackpos;
1042 term2 = stack[stackpos-1];
1043 --stackpos;
1044
1045 assert(func != fncvpower || SCIPexprGetOperator(term2) == SCIP_EXPR_CONST);
1046 assert(func != fnvcpower || SCIPexprGetOperator(term1) == SCIP_EXPR_CONST);
1047
1048 if( SCIPexprGetOperator(term1) == SCIP_EXPR_CONST )
1049 {
1050 /* use intpower if exponent is an integer constant, otherwise use realpower */
1051 if( SCIPisIntegral(scip, SCIPexprGetOpReal(term1)) )
1052 {
1053 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_INTPOWER, term2, (int)SCIPexprGetOpReal(term1)) );
1054 SCIPexprFreeDeep(blkmem, &term1);
1055 }
1056 else
1057 {
1058 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_REALPOWER, term2, SCIPexprGetOpReal(term1)) );
1059 SCIPexprFreeDeep(blkmem, &term1);
1060 }
1061 }
1062 else
1063 {
1064 /* term2^term1 = exp(log(term2)*term1) */
1065 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_LOG, term2) );
1066 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_MUL, e, term1) );
1067 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_EXP, e) );
1068 }
1069
1070 break;
1071 }
1072
1073 case fnsignpower: /* sign(x)*abs(x)^c */
1074 {
1075 SCIPdebugPrintf("signpower\n");
1076
1077 assert(stackpos >= 2);
1078 term1 = stack[stackpos-1];
1079 --stackpos;
1080 term2 = stack[stackpos-1];
1081 --stackpos;
1082
1083 if( SCIPexprGetOperator(term1) == SCIP_EXPR_CONST )
1084 {
1085 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_SIGNPOWER, term2, SCIPexprGetOpReal(term1)) );
1086 SCIPexprFreeDeep(blkmem, &term1);
1087 }
1088 else
1089 {
1090 SCIPerrorMessage("signpower with non-constant exponent not supported.\n");
1091 rc = SCIP_READERROR;
1092 goto TERMINATE;
1093 }
1094
1095 break;
1096 }
1097
1098 case fnpi:
1099 {
1100 SCIPdebugPrintf("pi\n");
1101
1102 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_CONST, M_PI) );
1103 break;
1104 }
1105
1106 case fndiv:
1107 {
1108 SCIPdebugPrintf("divide\n");
1109
1110 assert(stackpos >= 2);
1111 term1 = stack[stackpos-1];
1112 --stackpos;
1113 term2 = stack[stackpos-1];
1114 --stackpos;
1115
1116 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_DIV, term2, term1) );
1117 break;
1118 }
1119
1120 case fnabs:
1121 {
1122 SCIPdebugPrintf("abs\n");
1123
1124 assert(stackpos >= 1);
1125 term1 = stack[stackpos-1];
1126 --stackpos;
1127
1128 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_ABS, term1) );
1129 break;
1130 }
1131
1132 case fnpoly: /* univariate polynomial */
1133 {
1134 SCIPdebugPrintf("univariate polynomial of degree %d\n", nargs-2);
1135 assert(nargs >= 0);
1136 switch( nargs )
1137 {
1138 case 0:
1139 {
1140 term1 = stack[stackpos-1];
1141 --stackpos;
1142
1143 SCIPexprFreeDeep(blkmem, &term1);
1144 SCIP_CALL( SCIPexprCreate(blkmem, &e, SCIP_EXPR_CONST, 0.0) );
1145 break;
1146 }
1147
1148 case 1: /* "constant" polynomial */
1149 {
1150 e = stack[stackpos-1];
1151 --stackpos;
1152
1153 /* delete variable of polynomial */
1154 SCIPexprFreeDeep(blkmem, &stack[stackpos-1]);
1155 --stackpos;
1156
1157 break;
1158 }
1159
1160 default: /* polynomial is at least linear */
1161 {
1162 SCIP_EXPRDATA_MONOMIAL** monomials;
1163 SCIP_Real exponent;
1164 SCIP_Real constant;
1165 int nmonomials;
1166 int zero;
1167
1168 nmonomials = nargs-2;
1169 SCIP_CALL( SCIPallocBufferArray(scip, &monomials, nargs-2) );
1170
1171 zero = 0;
1172 constant = 0.0;
1173 for( ; nargs > 1; --nargs )
1174 {
1175 assert(stackpos > 0);
1176
1177 term1 = stack[stackpos-1];
1178 assert(SCIPexprGetOperator(term1) == SCIP_EXPR_CONST);
1179
1180 if( nargs > 2 )
1181 {
1182 exponent = (SCIP_Real)(nargs-2);
1183 SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[nargs-3], SCIPexprGetOpReal(term1), 1, &zero, &exponent) );
1184 }
1185 else
1186 constant = SCIPexprGetOpReal(term1);
1187
1188 SCIPexprFreeDeep(blkmem, &term1);
1189 --stackpos;
1190 }
1191
1192 assert(stackpos > 0);
1193 term1 = stack[stackpos-1];
1194 --stackpos;
1195
1196 SCIP_CALL( SCIPexprCreatePolynomial(blkmem, &e, 1, &term1, nmonomials, monomials, constant, FALSE) );
1197 SCIPfreeBufferArray(scip, &monomials);
1198 }
1199 }
1200 nargs = -1;
1201 break;
1202 }
1203
1204 /* @todo some of these we could also support */
1205 case fnerrf:
1206 case fnceil: case fnfloor: case fnround:
1207 case fnmod: case fntrunc: case fnsign:
1208 case fnarctan: case fndunfm:
1209 case fndnorm: case fnerror: case fnfrac: case fnerrorl:
1210 case fnfact /* factorial */:
1211 case fnunfmi /* uniform random number */:
1212 case fnncpf /* fischer: sqrt(x1^2+x2^2+2*x3) */:
1213 case fnncpcm /* chen-mangasarian: x1-x3*ln(1+exp((x1-x2)/x3))*/:
1214 case fnentropy /* x*ln(x) */: case fnsigmoid /* 1/(1+exp(-x)) */:
1215 case fnboolnot: case fnbooland:
1216 case fnboolor: case fnboolxor: case fnboolimp:
1217 case fnbooleqv: case fnrelopeq: case fnrelopgt:
1218 case fnrelopge: case fnreloplt: case fnrelople:
1219 case fnrelopne: case fnifthen:
1220 case fnedist /* euclidian distance */:
1221 case fncentropy /* x*ln((x+d)/(y+d))*/:
1222 case fngamma: case fnloggamma: case fnbeta:
1223 case fnlogbeta: case fngammareg: case fnbetareg:
1224 case fnsinh: case fncosh: case fntanh:
1225 case fnncpvusin /* veelken-ulbrich */:
1226 case fnncpvupow /* veelken-ulbrich */:
1227 case fnbinomial:
1228 case fnarccos:
1229 case fnarcsin: case fnarctan2 /* arctan(x2/x1) */:
1230 default :
1231 {
1232 SCIPdebugPrintf("nr. %d - unsupported. Error.\n", (int)func);
1233 SCIPinfoMessage(scip, NULL, "Error: GAMS function %s not supported.\n", GamsFuncCodeName[func]);
1234 rc = SCIP_READERROR;
1235 goto TERMINATE;
1236 }
1237 } /*lint !e788*/
1238 break;
1239 }
1240
1241 case nlEnd: /* end of instruction list */
1242 default:
1243 {
1244 SCIPinfoMessage(scip, NULL, "Error: GAMS opcode %s not supported.\n", GamsOpCodeName[opcode]);
1245 rc = SCIP_READERROR;
1246 goto TERMINATE;
1247 }
1248 } /*lint !e788*/
1249
1250 if( e != NULL )
1251 {
1252 if( stackpos >= stacksize )
1253 {
1254 stacksize = SCIPcalcMemGrowSize(scip, stackpos+1);
1255 SCIP_CALL( SCIPreallocBufferArray(scip, &stack, stacksize) );
1256 }
1257 assert(stackpos < stacksize);
1258 stack[stackpos] = e;
1259 ++stackpos;
1260 }
1261 }
1262
1263 /* there should be exactly one element on the stack, which will be the root of our expression tree */
1264 assert(stackpos == 1);
1265
1266 SCIP_CALL( SCIPexprtreeCreate(blkmem, exprtree, stack[0], nvars, 0, NULL) );
1267 SCIP_CALL( SCIPexprtreeSetVars(*exprtree, nvars, vars) );
1268
1269 TERMINATE:
1270 SCIPfreeBufferArray(scip, &vars);
1271 SCIPfreeBufferArray(scip, &stack);
1272 SCIPhashmapFree(&var2idx);
1273
1274 return rc;
1275 }
1276
1277 /** creates a SCIP problem from a GMO */
SCIPcreateProblemReaderGmo(SCIP * scip,gmoRec_t * gmo,const char * indicatorfile,int mipstart)1278 SCIP_RETCODE SCIPcreateProblemReaderGmo(
1279 SCIP* scip, /**< SCIP data structure */
1280 gmoRec_t* gmo, /**< GAMS Model Object */
1281 const char* indicatorfile, /**< name of file with indicator specification, or NULL */
1282 int mipstart /**< how to pass initial variable levels from GMO to SCIP */
1283 )
1284 {
1285 char buffer[GMS_SSSIZE];
1286 gevHandle_t gev;
1287 SCIP_Bool objnonlinear;
1288 SCIP_VAR** vars;
1289 SCIP_Real minprior;
1290 SCIP_Real maxprior;
1291 int i;
1292 SCIP_Real* coefs = NULL;
1293 int* indices = NULL;
1294 int* nlflag;
1295 SCIP_VAR** consvars = NULL;
1296 SCIP_VAR** quadvars1;
1297 SCIP_VAR** quadvars2;
1298 SCIP_Real* quadcoefs;
1299 int* qrow;
1300 int* qcol;
1301 SCIP_CONS* con;
1302 int numSos1, numSos2, nzSos;
1303 SCIP_PROBDATA* probdata;
1304 int* opcodes;
1305 int* fields;
1306 SCIP_Real* constants;
1307 int nindics;
1308 int* indicrows;
1309 int* indiccols;
1310 int* indiconvals;
1311 int indicidx;
1312 size_t namemem;
1313 SCIP_RETCODE rc = SCIP_OKAY;
1314
1315 assert(scip != NULL);
1316 assert(gmo != NULL);
1317
1318 gev = (gevHandle_t) gmoEnvironment(gmo);
1319 assert(gev != NULL);
1320
1321 /* we want a real objective function, if it is linear, otherwise keep the GAMS single-variable-objective? */
1322 gmoObjReformSet(gmo, 1);
1323 gmoObjStyleSet(gmo, (int)gmoObjType_Fun);
1324 #if 0
1325 if( gmoObjNLNZ(gmo) > 0 )
1326 gmoObjStyleSet(gmo, gmoObjType_Var);
1327 objnonlinear = FALSE;
1328 #else
1329 objnonlinear = gmoObjStyle(gmo) == (int)gmoObjType_Fun && gmoObjNLNZ(gmo) > 0;
1330 #endif
1331
1332 /* we want to start indexing at 0 */
1333 gmoIndexBaseSet(gmo, 0);
1334
1335 /* we want GMO to use SCIP's value for infinity */
1336 gmoPinfSet(gmo, SCIPinfinity(scip));
1337 gmoMinfSet(gmo, -SCIPinfinity(scip));
1338
1339 /* create SCIP problem */
1340 SCIP_CALL( SCIPallocMemory(scip, &probdata) );
1341 BMSclearMemory(probdata);
1342
1343 (void) gmoNameInput(gmo, buffer);
1344 SCIP_CALL( SCIPcreateProb(scip, buffer,
1345 probdataDelOrigGmo, probdataTransGmo, probdataDelTransGmo,
1346 probdataInitSolGmo, probdataExitSolGmo, probdataCopyGmo,
1347 probdata) );
1348
1349 /* initialize QMaker, if nonlinear */
1350 if( gmoNLNZ(gmo) > 0 || objnonlinear )
1351 gmoUseQSet(gmo, 1);
1352
1353 /* get data on indicator constraints from options object */
1354 nindics = 0;
1355 indicrows = NULL;
1356 indiccols = NULL;
1357 indiconvals = NULL;
1358 #ifndef WITH_GAMS
1359 if( indicatorfile != NULL && *indicatorfile != '\0' )
1360 {
1361 optHandle_t opt;
1362 int itype;
1363
1364 if( !optCreate(&opt, buffer, sizeof(buffer)) )
1365 {
1366 SCIPerrorMessage("*** Could not create optionfile handle: %s\n", buffer);
1367 return SCIP_ERROR;
1368 }
1369
1370 #if GMOAPIVERSION < 13
1371 (void) gevGetStrOpt(gev, gevNameSysDir, buffer);
1372 strcat(buffer, "optscip.def");
1373 if( optReadDefinition(opt, buffer) )
1374 #else
1375 if( optReadDefinitionFromPChar(opt, (char*)"indic indicator\ngeneral group 1 1 Dot options and indicators") )
1376 #endif
1377 {
1378 for( i = 1; i <= optMessageCount(opt); ++i )
1379 {
1380 optGetMessage(opt, i, buffer, &itype);
1381 if( itype <= (int) optMsgFileLeave || itype == (int) optMsgUserError )
1382 gevLogStat(gev, buffer);
1383 }
1384 optClearMessages(opt);
1385 return SCIP_ERROR;
1386 }
1387
1388 (void) optReadParameterFile(opt, indicatorfile);
1389 for( i = 1; i <= optMessageCount(opt); ++i )
1390 {
1391 optGetMessage(opt, i, buffer, &itype);
1392 if( itype <= (int) optMsgFileLeave || itype == (int) optMsgUserError )
1393 gevLogStat(gev, buffer);
1394 }
1395 optClearMessages(opt);
1396
1397 if( optIndicatorCount(opt, &i) > 0 )
1398 {
1399 SCIP_CALL( SCIPallocBufferArray(scip, &indicrows, gmoM(gmo)) );
1400 SCIP_CALL( SCIPallocBufferArray(scip, &indiccols, gmoM(gmo)) );
1401 SCIP_CALL( SCIPallocBufferArray(scip, &indiconvals, gmoM(gmo)) );
1402 if( gmoGetIndicatorMap(gmo, opt, 1, &nindics, indicrows, indiccols, indiconvals) != 0 )
1403 {
1404 SCIPerrorMessage("failed to get indicator mapping\n");
1405 return SCIP_ERROR;
1406 }
1407 }
1408
1409 (void) optFree(&opt);
1410 }
1411 #endif
1412 assert(indicrows != NULL || nindics == 0);
1413 assert(indiccols != NULL || nindics == 0);
1414 assert(indiconvals != NULL || nindics == 0);
1415
1416 namemem = (gmoN(gmo) + gmoM(gmo)) * sizeof(char*);
1417
1418 probdata->nvars = gmoN(gmo);
1419 SCIP_CALL( SCIPallocMemoryArray(scip, &probdata->vars, probdata->nvars) ); /*lint !e666*/
1420 vars = probdata->vars;
1421
1422 /* compute range of variable priorities */
1423 minprior = SCIPinfinity(scip);
1424 maxprior = 0.0;
1425 if( gmoPriorOpt(gmo) && gmoNDisc(gmo) > 0 )
1426 {
1427 for (i = 0; i < gmoN(gmo); ++i)
1428 {
1429 if( gmoGetVarTypeOne(gmo, i) == (int) gmovar_X )
1430 continue; /* GAMS forbids branching priorities for continuous variables */
1431 if( gmoGetVarPriorOne(gmo,i) < minprior )
1432 minprior = gmoGetVarPriorOne(gmo,i);
1433 if( gmoGetVarPriorOne(gmo,i) > maxprior )
1434 maxprior = gmoGetVarPriorOne(gmo,i);
1435 }
1436 }
1437
1438 /* get objective functions coefficients */
1439 SCIP_CALL( SCIPallocBufferArray(scip, &coefs, gmoN(gmo)+1) ); /* +1 if we have to transform the objective into a constraint */
1440 if( !objnonlinear )
1441 {
1442 if( gmoObjStyle(gmo) == (int) gmoObjType_Fun )
1443 (void) gmoGetObjVector(gmo, coefs, NULL);
1444 else
1445 coefs[gmoObjVar(gmo)] = 1.0;
1446 }
1447 else
1448 {
1449 BMSclearMemoryArray(coefs, gmoN(gmo));
1450 }
1451
1452 /* add variables */
1453 for( i = 0; i < gmoN(gmo); ++i )
1454 {
1455 SCIP_VARTYPE vartype;
1456 SCIP_Real lb, ub;
1457 lb = gmoGetVarLowerOne(gmo, i);
1458 ub = gmoGetVarUpperOne(gmo, i);
1459 switch( gmoGetVarTypeOne(gmo, i) )
1460 {
1461 case gmovar_SC:
1462 lb = 0.0;
1463 /*lint -fallthrough*/
1464 case gmovar_X:
1465 case gmovar_S1:
1466 case gmovar_S2:
1467 vartype = SCIP_VARTYPE_CONTINUOUS;
1468 break;
1469 case gmovar_B:
1470 vartype = SCIP_VARTYPE_BINARY;
1471 break;
1472 case gmovar_SI:
1473 lb = 0.0;
1474 /*lint -fallthrough*/
1475 case gmovar_I:
1476 vartype = SCIP_VARTYPE_INTEGER;
1477 break;
1478 default:
1479 SCIPerrorMessage("Unknown variable type.\n");
1480 return SCIP_INVALIDDATA;
1481 }
1482 if( gmoDict(gmo) )
1483 {
1484 (void) gmoGetVarNameOne(gmo, i, buffer);
1485 if( nindics == 0 )
1486 namemem += strlen(buffer) + 1;
1487 }
1488 else
1489 sprintf(buffer, "x%d", i);
1490 SCIP_CALL( SCIPcreateVar(scip, &vars[i], buffer, lb, ub, coefs[i], vartype, TRUE, FALSE, NULL, NULL, NULL, NULL, NULL) );
1491 SCIP_CALL( SCIPaddVar(scip, vars[i]) );
1492 SCIPdebugMessage("added variable ");
1493 SCIPdebug( SCIPprintVar(scip, vars[i], NULL) );
1494
1495 if( gmoPriorOpt(gmo) && minprior < maxprior && gmoGetVarTypeOne(gmo, i) != (int) gmovar_X )
1496 {
1497 /* in GAMS: higher priorities are given by smaller .prior values
1498 in SCIP: variables with higher branch priority are always preferred to variables with lower priority in selection of branching variable
1499 thus, we scale the values from GAMS to lie between 0 (lowest prior) and 1000 (highest prior)
1500 */
1501 int branchpriority = (int)(1000.0 / (maxprior - minprior) * (maxprior - gmoGetVarPriorOne(gmo, i)));
1502 SCIP_CALL( SCIPchgVarBranchPriority(scip, vars[i], branchpriority) );
1503 }
1504 }
1505
1506 /* setup bound disjunction constraints for semicontinuous/semiinteger variables by saying x <= 0 or x >= gmoGetVarLower */
1507 if( gmoGetVarTypeCnt(gmo, (int) gmovar_SC) || gmoGetVarTypeCnt(gmo, (int) gmovar_SI) )
1508 {
1509 SCIP_BOUNDTYPE bndtypes[2];
1510 SCIP_Real bnds[2];
1511 SCIP_VAR* bndvars[2];
1512 SCIP_CONS* cons;
1513 char name[SCIP_MAXSTRLEN];
1514
1515 bndtypes[0] = SCIP_BOUNDTYPE_UPPER;
1516 bndtypes[1] = SCIP_BOUNDTYPE_LOWER;
1517 bnds[0] = 0.0;
1518
1519 for( i = 0; i < gmoN(gmo); ++i )
1520 {
1521 if( gmoGetVarTypeOne(gmo, i) != (int) gmovar_SC && gmoGetVarTypeOne(gmo, i) != (int) gmovar_SI )
1522 continue;
1523
1524 bnds[1] = gmoGetVarLowerOne(gmo, i);
1525 /* skip bound disjunction if lower bound is 0 (if continuous) or 1 (if integer)
1526 * since this is a trivial disjunction, but would raise an assert in SCIPcreateConsBounddisjunction
1527 */
1528 if( gmoGetVarTypeOne(gmo, i) == (int) gmovar_SC && !SCIPisPositive(scip, bnds[1]) )
1529 continue;
1530 if( gmoGetVarTypeOne(gmo, i) == (int) gmovar_SI && !SCIPisPositive(scip, bnds[1]-1.0) )
1531 continue;
1532
1533 bndvars[0] = vars[i];
1534 bndvars[1] = vars[i];
1535 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "semi%s_%s", gmoGetVarTypeOne(gmo, i) == (int) gmovar_SC ? "con" : "int", SCIPvarGetName(vars[i]));
1536 SCIP_CALL( SCIPcreateConsBounddisjunction(scip, &cons, name, 2, bndvars, bndtypes, bnds,
1537 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1538 SCIP_CALL( SCIPaddCons(scip, cons) );
1539 SCIPdebugMessage("added constraint ");
1540 SCIPdebug( SCIPprintCons(scip, cons, NULL) );
1541 SCIP_CALL( SCIPreleaseCons(scip, &cons) );
1542 }
1543 }
1544
1545 SCIP_CALL( SCIPallocBufferArray(scip, &consvars, gmoN(gmo)+1) ); /* +1 if we have to transform the objective into a constraint */
1546
1547 /* setup SOS constraints */
1548 gmoGetSosCounts(gmo, &numSos1, &numSos2, &nzSos);
1549 if( nzSos > 0 )
1550 {
1551 int numSos;
1552 int* sostype;
1553 int* sosbeg;
1554 int* sosind;
1555 double* soswt;
1556 int j, k;
1557
1558 numSos = numSos1 + numSos2;
1559 SCIP_CALL( SCIPallocBufferArray(scip, &sostype, numSos) );
1560 SCIP_CALL( SCIPallocBufferArray(scip, &sosbeg, numSos+1) );
1561 SCIP_CALL( SCIPallocBufferArray(scip, &sosind, nzSos) );
1562 SCIP_CALL( SCIPallocBufferArray(scip, &soswt, nzSos) );
1563
1564 (void) gmoGetSosConstraints(gmo, sostype, sosbeg, sosind, soswt);
1565
1566 for( i = 0; i < numSos; ++i )
1567 {
1568 for( j = sosbeg[i], k = 0; j < sosbeg[i+1]; ++j, ++k )
1569 {
1570 consvars[k] = vars[sosind[j]];
1571 assert(gmoGetVarTypeOne(gmo, sosind[j]) == (sostype[i] == 1 ? (int) gmovar_S1 : (int) gmovar_S2));
1572 }
1573
1574 sprintf(buffer, "sos%d", i);
1575 if( sostype[i] == 1 )
1576 {
1577 SCIP_CALL( SCIPcreateConsSOS1(scip, &con, buffer, k, consvars, &soswt[sosbeg[i]], TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
1578 }
1579 else
1580 {
1581 assert(sostype[i] == 2);
1582 SCIP_CALL( SCIPcreateConsSOS2(scip, &con, buffer, k, consvars, &soswt[sosbeg[i]], TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
1583 }
1584
1585 SCIP_CALL( SCIPaddCons(scip, con) );
1586 SCIPdebugMessage("added constraint ");
1587 SCIPdebug( SCIPprintCons(scip, con, NULL) );
1588 SCIP_CALL( SCIPreleaseCons(scip, &con) );
1589 }
1590
1591 SCIPfreeBufferArray(scip, &sostype);
1592 SCIPfreeBufferArray(scip, &sosbeg);
1593 SCIPfreeBufferArray(scip, &sosind);
1594 SCIPfreeBufferArray(scip, &soswt);
1595 }
1596
1597 /* setup regular constraints */
1598 SCIP_CALL( SCIPallocBufferArray(scip, &indices, gmoN(gmo)) );
1599 indicidx = 0;
1600
1601 /* alloc some memory, if nonlinear */
1602 if( gmoNLNZ(gmo) > 0 || objnonlinear )
1603 {
1604 SCIP_CALL( SCIPallocBufferArray(scip, &nlflag, gmoN(gmo)) );
1605
1606 SCIP_CALL( SCIPallocBufferArray(scip, &quadvars1, gmoMaxQNZ(gmo)) );
1607 SCIP_CALL( SCIPallocBufferArray(scip, &quadvars2, gmoMaxQNZ(gmo)) );
1608 SCIP_CALL( SCIPallocBufferArray(scip, &quadcoefs, gmoMaxQNZ(gmo)) );
1609 SCIP_CALL( SCIPallocBufferArray(scip, &qrow, gmoMaxQNZ(gmo)) );
1610 SCIP_CALL( SCIPallocBufferArray(scip, &qcol, gmoMaxQNZ(gmo)) );
1611
1612 SCIP_CALL( SCIPallocBufferArray(scip, &opcodes, gmoNLCodeSizeMaxRow(gmo)+1) );
1613 SCIP_CALL( SCIPallocBufferArray(scip, &fields, gmoNLCodeSizeMaxRow(gmo)+1) );
1614 SCIP_CALL( SCIPduplicateBufferArray(scip, &constants, (double*)gmoPPool(gmo), gmoNLConst(gmo)) );
1615
1616 /* translate special GAMS constants into SCIP variants (gmo does not seem to do this...) */
1617 for( i = 0; i < gmoNLConst(gmo); ++i )
1618 {
1619 if( constants[i] == GMS_SV_PINF )
1620 constants[i] = SCIPinfinity(scip);
1621 else if( constants[i] == GMS_SV_MINF )
1622 constants[i] = -SCIPinfinity(scip);
1623 else if( constants[i] == GMS_SV_EPS )
1624 constants[i] = 0.0;
1625 else if( constants[i] == GMS_SV_UNDEF || constants[i] == GMS_SV_NA || constants[i] == GMS_SV_NAINT || constants[i] == GMS_SV_ACR )
1626 {
1627 SCIPwarningMessage(scip, "Constant %e in nonlinear expressions constants pool cannot be handled by SCIP.\n", constants[i]);
1628 constants[i] = SCIP_INVALID;
1629 }
1630 else if( constants[i] <= -SCIPinfinity(scip) )
1631 constants[i] = -SCIPinfinity(scip);
1632 else if( constants[i] >= SCIPinfinity(scip) )
1633 constants[i] = SCIPinfinity(scip);
1634 }
1635 }
1636 else
1637 {
1638 nlflag = NULL;
1639
1640 quadvars1 = NULL;
1641 quadvars2 = NULL;
1642 quadcoefs = NULL;
1643 qrow = NULL;
1644 qcol = NULL;
1645
1646 opcodes = NULL;
1647 fields = NULL;
1648 constants = NULL;
1649 }
1650
1651 for( i = 0; i < gmoM(gmo); ++i )
1652 {
1653 double lhs;
1654 double rhs;
1655 switch( gmoGetEquTypeOne(gmo, i) )
1656 {
1657 case gmoequ_E:
1658 lhs = rhs = gmoGetRhsOne(gmo, i);
1659 break;
1660 case gmoequ_G:
1661 lhs = gmoGetRhsOne(gmo, i);
1662 rhs = SCIPinfinity(scip);
1663 break;
1664 case gmoequ_L:
1665 lhs = -SCIPinfinity(scip);
1666 rhs = gmoGetRhsOne(gmo, i);
1667 break;
1668 case gmoequ_N:
1669 lhs = -SCIPinfinity(scip);
1670 rhs = SCIPinfinity(scip);
1671 break;
1672 case gmoequ_X:
1673 SCIPerrorMessage("External functions not supported by SCIP.\n");
1674 return SCIP_INVALIDDATA;
1675 case gmoequ_C:
1676 SCIPerrorMessage("Conic constraints not supported by SCIP interface yet.\n");
1677 return SCIP_INVALIDDATA;
1678 case gmoequ_B:
1679 SCIPerrorMessage("Logic constraints not supported by SCIP interface yet.\n");
1680 return SCIP_INVALIDDATA;
1681 default:
1682 SCIPerrorMessage("unknown equation type.\n");
1683 return SCIP_INVALIDDATA;
1684 }
1685
1686 if( gmoDict(gmo) )
1687 {
1688 (void) gmoGetEquNameOne(gmo, i, buffer);
1689 if( nindics == 0 )
1690 namemem += strlen(buffer) + 1;
1691 }
1692 else
1693 sprintf(buffer, "e%d", i);
1694
1695 con = NULL;
1696 switch( gmoGetEquOrderOne(gmo, i) )
1697 {
1698 case gmoorder_L:
1699 {
1700 /* linear constraint */
1701 int j, nz, nlnz;
1702 (void) gmoGetRowSparse(gmo, i, indices, coefs, NULL, &nz, &nlnz);
1703 assert(nlnz == 0);
1704
1705 for( j = 0; j < nz; ++j )
1706 consvars[j] = vars[indices[j]];
1707
1708 /* create indicator constraint, if we are at one */
1709 if( indicidx < nindics && indicrows[indicidx] == i ) /*lint !e613*/
1710 {
1711 SCIP_VAR* binvar;
1712
1713 binvar = vars[indiccols[indicidx]]; /*lint !e613*/
1714 if( SCIPvarGetType(binvar) != SCIP_VARTYPE_BINARY )
1715 {
1716 SCIPerrorMessage("Indicator variable <%s> is not of binary type.\n", SCIPvarGetName(binvar));
1717 return SCIP_ERROR;
1718 }
1719
1720 assert(indiconvals[indicidx] == 0 || indiconvals[indicidx] == 1); /*lint !e613*/
1721 if( indiconvals[indicidx] == 0 ) /*lint !e613*/
1722 {
1723 SCIP_CALL( SCIPgetNegatedVar(scip, binvar, &binvar) );
1724 }
1725
1726 if( !SCIPisInfinity(scip, rhs) )
1727 {
1728 SCIP_CALL( SCIPcreateConsIndicator(scip, &con, buffer, binvar, nz, consvars, coefs, rhs,
1729 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
1730
1731 if( !SCIPisInfinity(scip, -lhs) )
1732 {
1733 SCIP_CALL( SCIPaddCons(scip, con) );
1734 SCIPdebugMessage("added constraint ");
1735 SCIPdebug( SCIPprintCons(scip, con, NULL) );
1736 SCIP_CALL( SCIPreleaseCons(scip, &con) );
1737 con = NULL;
1738 }
1739 }
1740 if( !SCIPisInfinity(scip, -lhs) )
1741 {
1742 for( j = 0; j < nz; ++j )
1743 coefs[j] = -coefs[j];
1744 SCIP_CALL( SCIPcreateConsIndicator(scip, &con, buffer, binvar, nz, consvars, coefs, -lhs,
1745 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
1746 }
1747
1748 ++indicidx;
1749 }
1750 else
1751 {
1752 SCIP_CALL( SCIPcreateConsLinear(scip, &con, buffer, nz, consvars, coefs, lhs, rhs,
1753 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1754 }
1755
1756 break;
1757 }
1758
1759 case gmoorder_Q:
1760 {
1761 /* quadratic constraint */
1762 int j, nz, nlnz, qnz;
1763
1764 assert(qcol != NULL);
1765 assert(qrow != NULL);
1766 assert(quadcoefs != NULL);
1767 assert(quadvars1 != NULL);
1768 assert(quadvars2 != NULL);
1769
1770 (void) gmoGetRowSparse(gmo, i, indices, coefs, NULL, &nz, &nlnz);
1771 for( j = 0; j < nz; ++j )
1772 consvars[j] = vars[indices[j]];
1773
1774 qnz = gmoGetRowQNZOne(gmo,i);
1775 (void) gmoGetRowQ(gmo, i, qcol, qrow, quadcoefs);
1776 for( j = 0; j < qnz; ++j )
1777 {
1778 assert(qcol[j] >= 0);
1779 assert(qrow[j] >= 0);
1780 assert(qcol[j] < gmoN(gmo));
1781 assert(qrow[j] < gmoN(gmo));
1782 quadvars1[j] = vars[qcol[j]];
1783 quadvars2[j] = vars[qrow[j]];
1784 if( qcol[j] == qrow[j] )
1785 quadcoefs[j] /= 2.0; /* for some strange reason, the coefficients on the diagonal are multiplied by 2 in GMO. */
1786 }
1787
1788 SCIP_CALL( SCIPcreateConsQuadratic(scip, &con, buffer, nz, consvars, coefs, qnz, quadvars1, quadvars2, quadcoefs, lhs, rhs,
1789 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
1790 break;
1791 }
1792
1793 case gmoorder_NL:
1794 {
1795 /* nonlinear constraint */
1796 int j, nz, nlnz, linnz;
1797 int codelen;
1798 SCIP_EXPRTREE* exprtree;
1799
1800 assert(nlflag != NULL);
1801
1802 (void) gmoGetRowSparse(gmo, i, indices, coefs, nlflag, &nz, &nlnz);
1803 linnz = 0;
1804 for( j = 0; j < nz; ++j )
1805 {
1806 if( !nlflag[j] )
1807 {
1808 consvars[linnz] = vars[indices[j]];
1809 coefs[linnz] = coefs[j];
1810 ++linnz;
1811 }
1812 }
1813
1814 (void) gmoDirtyGetRowFNLInstr(gmo, i, &codelen, opcodes, fields);
1815 rc = makeExprtree(scip, gmo, codelen, opcodes, fields, constants, &exprtree);
1816 if( rc == SCIP_READERROR )
1817 {
1818 SCIPinfoMessage(scip, NULL, "Error processing nonlinear instructions of equation %s.\n", buffer);
1819 goto TERMINATE;
1820 }
1821 SCIP_CALL( rc );
1822
1823 SCIP_CALL( SCIPcreateConsNonlinear(scip, &con, buffer, linnz, consvars, coefs, 1, &exprtree, NULL, lhs, rhs,
1824 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1825
1826 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
1827 break;
1828 }
1829
1830 default:
1831 SCIPerrorMessage("Unexpected equation order.\n");
1832 return SCIP_INVALIDDATA;
1833 }
1834
1835 assert(con != NULL);
1836 SCIP_CALL( SCIPaddCons(scip, con) );
1837 SCIPdebugMessage("added constraint ");
1838 SCIPdebug( SCIPprintCons(scip, con, NULL) );
1839 SCIP_CALL( SCIPreleaseCons(scip, &con) );
1840
1841 /* @todo do something about this */
1842 if( indicidx < nindics && indicrows[indicidx] == i ) /*lint !e613*/
1843 {
1844 SCIPerrorMessage("Only linear constraints can be indicatored, currently.\n");
1845 return SCIP_ERROR;
1846 }
1847 }
1848
1849 if( objnonlinear )
1850 {
1851 /* make constraint out of nonlinear objective function */
1852 int j, nz, nlnz, qnz;
1853 double lhs, rhs;
1854
1855 assert(gmoGetObjOrder(gmo) == (int) gmoorder_L || gmoGetObjOrder(gmo) == (int) gmoorder_Q || gmoGetObjOrder(gmo) == (int) gmoorder_NL);
1856
1857 SCIP_CALL( SCIPcreateVar(scip, &probdata->objvar, "xobj", -SCIPinfinity(scip), SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE, NULL, NULL, NULL, NULL, NULL) );
1858 SCIP_CALL( SCIPaddVar(scip, probdata->objvar) );
1859 SCIPdebugMessage("added objective variable ");
1860 SCIPdebug( SCIPprintVar(scip, probdata->objvar, NULL) );
1861
1862 if( gmoGetObjOrder(gmo) != (int) gmoorder_NL )
1863 {
1864 assert(qcol != NULL);
1865 assert(qrow != NULL);
1866 assert(quadcoefs != NULL);
1867 assert(quadvars1 != NULL);
1868 assert(quadvars2 != NULL);
1869
1870 (void) gmoGetObjSparse(gmo, indices, coefs, NULL, &nz, &nlnz);
1871 for( j = 0; j < nz; ++j )
1872 consvars[j] = vars[indices[j]];
1873
1874 consvars[nz] = probdata->objvar;
1875 coefs[nz] = -1.0;
1876 ++nz;
1877
1878 qnz = gmoObjQNZ(gmo);
1879 (void) gmoGetObjQ(gmo, qcol, qrow, quadcoefs);
1880 for( j = 0; j < qnz; ++j )
1881 {
1882 assert(qcol[j] >= 0);
1883 assert(qrow[j] >= 0);
1884 assert(qcol[j] < gmoN(gmo));
1885 assert(qrow[j] < gmoN(gmo));
1886 quadvars1[j] = vars[qcol[j]];
1887 quadvars2[j] = vars[qrow[j]];
1888 if( qcol[j] == qrow[j] )
1889 quadcoefs[j] /= 2.0; /* for some strange reason, the coefficients on the diagonal are multiplied by 2 in GMO */
1890 }
1891
1892 if( gmoSense(gmo) == (int) gmoObj_Min )
1893 {
1894 lhs = -SCIPinfinity(scip);
1895 rhs = -gmoObjConst(gmo);
1896 }
1897 else
1898 {
1899 lhs = -gmoObjConst(gmo);
1900 rhs = SCIPinfinity(scip);
1901 }
1902
1903 SCIP_CALL( SCIPcreateConsQuadratic(scip, &con, "objective", nz, consvars, coefs, qnz, quadvars1, quadvars2, quadcoefs, lhs, rhs,
1904 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
1905 }
1906 else
1907 {
1908 SCIP_Real objfactor;
1909 int linnz;
1910 int codelen;
1911 SCIP_EXPRTREE* exprtree;
1912
1913 assert(nlflag != NULL);
1914
1915 (void) gmoGetObjSparse(gmo, indices, coefs, nlflag, &nz, &nlnz);
1916 linnz = 0;
1917 for( j = 0; j < nz; ++j )
1918 {
1919 if( !nlflag[j] )
1920 {
1921 coefs[linnz] = coefs[j];
1922 consvars[linnz] = vars[indices[j]];
1923 ++linnz;
1924 }
1925 }
1926
1927 consvars[linnz] = probdata->objvar;
1928 coefs[linnz] = -1.0;
1929 ++linnz;
1930
1931 objfactor = -1.0 / gmoObjJacVal(gmo);
1932
1933 (void) gmoDirtyGetObjFNLInstr(gmo, &codelen, opcodes, fields);
1934 rc = makeExprtree(scip, gmo, codelen, opcodes, fields, constants, &exprtree);
1935 if( rc == SCIP_READERROR )
1936 {
1937 SCIPinfoMessage(scip, NULL, "Error processing nonlinear instructions of objective %s.\n", gmoGetObjName(gmo, buffer));
1938 goto TERMINATE;
1939 }
1940 SCIP_CALL( rc );
1941
1942 if( gmoSense(gmo) == (int) gmoObj_Min )
1943 {
1944 lhs = -SCIPinfinity(scip);
1945 rhs = -gmoObjConst(gmo);
1946 }
1947 else
1948 {
1949 lhs = -gmoObjConst(gmo);
1950 rhs = SCIPinfinity(scip);
1951 }
1952
1953 SCIP_CALL( SCIPcreateConsNonlinear(scip, &con, "objective", linnz, consvars, coefs, 1, &exprtree, &objfactor, lhs, rhs,
1954 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1955
1956 SCIP_CALL( SCIPexprtreeFree(&exprtree) );
1957 }
1958
1959 SCIP_CALL( SCIPaddCons(scip, con) );
1960 SCIPdebugMessage("added objective constraint ");
1961 SCIPdebug( SCIPprintCons(scip, con, NULL) );
1962 SCIP_CALL( SCIPreleaseCons(scip, &con) );
1963 }
1964 else if( !SCIPisZero(scip, gmoObjConst(gmo)) )
1965 {
1966 /* handle constant term in linear objective by adding a fixed variable */
1967 SCIP_CALL( SCIPcreateVar(scip, &probdata->objconst, "objconst", 1.0, 1.0, gmoObjConst(gmo), SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE, NULL, NULL, NULL, NULL, NULL) );
1968 SCIP_CALL( SCIPaddVar(scip, probdata->objconst) );
1969 SCIPdebugMessage("added variable for objective constant: ");
1970 SCIPdebug( SCIPprintVar(scip, probdata->objconst, NULL) );
1971 }
1972
1973 if( gmoSense(gmo) == (int) gmoObj_Max )
1974 SCIP_CALL( SCIPsetObjsense(scip, SCIP_OBJSENSE_MAXIMIZE) );
1975
1976 /* set objective limit, if enabled */
1977 if( gevGetIntOpt(gev, gevUseCutOff) )
1978 {
1979 SCIP_CALL( SCIPsetObjlimit(scip, gevGetDblOpt(gev, gevCutOff)) );
1980 }
1981
1982 /* handle initial solution values */
1983 switch( mipstart )
1984 {
1985 case 0 :
1986 {
1987 /* don't pass any initial values to SCIP */
1988 break;
1989 }
1990
1991 case 2:
1992 case 3:
1993 {
1994 /* pass all initial values to SCIP and let SCIP check feasibility (2) or repair (3)
1995 * NOTE: mipstart=3 does not work as expected: heur_completesol does not run if values for all vars are given and for all integer variables integral values are given
1996 */
1997 SCIP_SOL* sol;
1998 SCIP_Real* vals;
1999 SCIP_Bool stored;
2000
2001 if( mipstart == 2 )
2002 {
2003 /* with this, SCIP will only check feasibility */
2004 SCIP_CALL( SCIPcreateOrigSol(scip, &sol, NULL) );
2005 }
2006 else
2007 {
2008 /* with this, SCIP will try to find a feasible solution close by to the initial values */
2009 SCIP_CALL( SCIPcreatePartialSol(scip, &sol, NULL) );
2010 }
2011
2012 SCIP_CALL( SCIPallocBufferArray(scip, &vals, gmoN(gmo)) );
2013 (void) gmoGetVarL(gmo, vals);
2014
2015 SCIP_CALL( SCIPsetSolVals(scip, sol, gmoN(gmo), probdata->vars, vals) );
2016
2017 /* if we have extra variable for objective, then need to set its value too */
2018 if( probdata->objvar != NULL )
2019 {
2020 double objval;
2021 int numErr;
2022 (void) gmoEvalFuncObj(gmo, vals, &objval, &numErr);
2023 if( numErr == 0 )
2024 {
2025 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->objvar, objval) );
2026 }
2027 }
2028
2029 /* if we have extra variable for objective constant, then need to set its value to 1.0 here too */
2030 if( probdata->objconst != NULL )
2031 {
2032 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->objconst, 1.0) );
2033 }
2034
2035 SCIP_CALL( SCIPaddSolFree(scip, &sol, &stored) );
2036 assert(stored);
2037
2038 SCIPfreeBufferArray(scip, &vals);
2039
2040 if( mipstart == 3 )
2041 {
2042 SCIPinfoMessage(scip, NULL, "Passed partial solution with values for all variables to SCIP.");
2043 }
2044
2045 break;
2046 }
2047
2048 case 1:
2049 case 4:
2050 {
2051 /* pass some initial value to SCIP and let SCIP complete solution */
2052 SCIP_SOL* sol;
2053 SCIP_Bool stored;
2054 double tryint = 0.0;
2055 int nknown;
2056
2057 if( mipstart == 4 )
2058 tryint = gevGetDblOpt(gev, gevTryInt);
2059
2060 SCIP_CALL( SCIPcreatePartialSol(scip, &sol, NULL) );
2061
2062 nknown = 0;
2063 for( i = 0; i < gmoN(gmo); ++i )
2064 {
2065 if( mipstart == 1 && (gmoGetVarTypeOne(gmo, i) == gmovar_B || gmoGetVarTypeOne(gmo, i) == gmovar_I || gmoGetVarTypeOne(gmo, i) == gmovar_SI) )
2066 {
2067 /* 1: set all integer variables */
2068 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->vars[i], gmoGetVarLOne(gmo, i)) );
2069 ++nknown;
2070 }
2071
2072 if( mipstart == 4 && (gmoGetVarTypeOne(gmo, i) == gmovar_B || gmoGetVarTypeOne(gmo, i) == gmovar_I || gmoGetVarTypeOne(gmo, i) == gmovar_SI) )
2073 {
2074 /* 4: set only integer variables with level close to an integral value, closeness decided by tryint */
2075 SCIP_Real val;
2076
2077 val = gmoGetVarLOne(gmo, i);
2078 if( fabs(round(val)-val) <= tryint )
2079 {
2080 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->vars[i], val) );
2081 ++nknown;
2082 }
2083 }
2084 }
2085
2086 /* if we have extra variable for objective constant, then can set its value to 1.0 here too */
2087 if( probdata->objconst != NULL )
2088 {
2089 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->objconst, 1.0) );
2090 }
2091
2092 SCIP_CALL( SCIPaddSolFree(scip, &sol, &stored) );
2093 assert(stored);
2094
2095 SCIPinfoMessage(scip, NULL, "Passed partial solution with values for %d variables (%.1f%%) to SCIP.", nknown, 100.0*(double)nknown/gmoN(gmo));
2096
2097 break;
2098 }
2099
2100 default:
2101 {
2102 SCIPwarningMessage(scip, "Setting mipstart = %d not supported. Ignored.\n", mipstart);
2103 }
2104 }
2105
2106 if( namemem > 1024 * 1024 && nindics == 0 )
2107 {
2108 namemem <<= 1; /* transformed problem has copy of names, so duplicate estimate */
2109 SCIPinfoMessage(scip, NULL, "Space for names approximately %0.2f MB. Use statement '<modelname>.dictfile=0;' to turn dictionary off.\n", namemem/(1024.0*1024.0));
2110 }
2111
2112 TERMINATE:
2113 SCIPfreeBufferArrayNull(scip, &coefs);
2114 SCIPfreeBufferArrayNull(scip, &indices);
2115 SCIPfreeBufferArrayNull(scip, &consvars);
2116 SCIPfreeBufferArrayNull(scip, &nlflag);
2117 SCIPfreeBufferArrayNull(scip, &quadvars1);
2118 SCIPfreeBufferArrayNull(scip, &quadvars2);
2119 SCIPfreeBufferArrayNull(scip, &quadcoefs);
2120 SCIPfreeBufferArrayNull(scip, &qrow);
2121 SCIPfreeBufferArrayNull(scip, &qcol);
2122 SCIPfreeBufferArrayNull(scip, &opcodes);
2123 SCIPfreeBufferArrayNull(scip, &fields);
2124 SCIPfreeBufferArrayNull(scip, &constants);
2125 SCIPfreeBufferArrayNull(scip, &indicrows);
2126 SCIPfreeBufferArrayNull(scip, &indiccols);
2127 SCIPfreeBufferArrayNull(scip, &indiconvals);
2128
2129 /* deinitialize QMaker, if nonlinear */
2130 if( gmoNLNZ(gmo) > 0 || objnonlinear )
2131 gmoUseQSet(gmo, 0);
2132
2133 return rc;
2134 }
2135
2136 /** check solution for feasibility and resolves by NLP solver, if necessary and possible */
2137 static
checkAndRepairSol(SCIP * scip,SCIP_Real * solvals,SCIP_Real * objval,SCIP_Bool resolvenlp,SCIP_Bool * success)2138 SCIP_RETCODE checkAndRepairSol(
2139 SCIP* scip, /**< SCIP data structure */
2140 SCIP_Real* solvals, /**< solution to check */
2141 SCIP_Real* objval, /**< objective value corresponding to solvals */
2142 SCIP_Bool resolvenlp, /**< whether NLP resolving is allowed */
2143 SCIP_Bool* success /**< to store whether solution is feasible or could be made feasible */
2144 )
2145 {
2146 SCIP_PROBDATA* probdata;
2147 SCIP_HEUR* heursubnlp;
2148 SCIP_SOL* sol;
2149
2150 assert(scip != NULL);
2151 assert(solvals != NULL);
2152 assert(success != NULL);
2153 assert(objval != NULL);
2154
2155 if( SCIPisTransformed(scip) )
2156 {
2157 /* cannot create solutions in SOLVED stage */
2158 SCIP_CALL( SCIPfreeTransform(scip) );
2159 }
2160
2161 probdata = SCIPgetProbData(scip);
2162 assert(probdata != NULL);
2163
2164 SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) );
2165 SCIP_CALL( SCIPsetSolVals(scip, sol, probdata->nvars, probdata->vars, solvals) );
2166 if( probdata->objvar != NULL )
2167 {
2168 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->objvar, *objval) );
2169 }
2170 if( probdata->objconst != NULL )
2171 {
2172 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->objconst, 1.0) );
2173 }
2174
2175 SCIP_CALL( SCIPcheckSolOrig(scip, sol, success, FALSE, FALSE) );
2176
2177 SCIP_CALL( SCIPfreeSol(scip, &sol) );
2178
2179 if( *success || !resolvenlp )
2180 return SCIP_OKAY;
2181
2182 /* assert that we checked already that resolving is possible and makes sense */
2183 assert(SCIPgetNContVars(scip) > 0);
2184 assert(SCIPgetNNlpis(scip) > 0);
2185
2186 heursubnlp = SCIPfindHeur(scip, "subnlp");
2187 assert(heursubnlp != NULL);
2188
2189 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "Attempt solving NLP from original problem with fixed discrete variables.\n");
2190
2191 /* create transformed problem and recreate sol in transformed problem, so subnlp heuristic can return result in it */
2192 SCIP_CALL( SCIPtransformProb(scip) );
2193 SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) );
2194 SCIP_CALL( SCIPsetSolVals(scip, sol, probdata->nvars, probdata->vars, solvals) );
2195 if( probdata->objvar != NULL )
2196 {
2197 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->objvar, *objval) );
2198 }
2199 if( probdata->objconst != NULL )
2200 {
2201 SCIP_CALL( SCIPsetSolVal(scip, sol, probdata->objconst, 1.0) );
2202 }
2203
2204 SCIP_CALL( SCIPresolveSolHeurSubNlp(scip, heursubnlp, sol, success, 100LL, 10.0) );
2205
2206 if( *success )
2207 {
2208 SCIP_CALL( SCIPcheckSolOrig(scip, sol, success, SCIPgetVerbLevel(scip) >= SCIP_VERBLEVEL_FULL, FALSE) );
2209
2210 if( *success )
2211 {
2212 SCIP_CALL( SCIPgetSolVals(scip, sol, probdata->nvars, probdata->vars, solvals) );
2213 *objval = SCIPgetSolOrigObj(scip, sol);
2214
2215 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "NLP solution is feasible, objective value = %.15e.\n", *objval);
2216 }
2217 else
2218 {
2219 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "NLP solution still not feasible, objective value = %.15e.\n", SCIPgetSolOrigObj(scip, sol));
2220 }
2221 }
2222 else
2223 {
2224 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "Failed to resolve NLP.\n");
2225 }
2226
2227 SCIP_CALL( SCIPfreeSol(scip, &sol) );
2228
2229 return SCIP_OKAY;
2230 }
2231
2232 /** stores solve information (solution, statistics) in a GMO */
2233 static
writeGmoSolution(SCIP * scip,gmoHandle_t gmo)2234 SCIP_RETCODE writeGmoSolution(
2235 SCIP* scip, /**< SCIP data structure */
2236 gmoHandle_t gmo /**< GAMS Model Object */
2237 )
2238 {
2239 SCIP_PROBDATA* probdata;
2240 int nrsol;
2241 SCIP_Real dualbound;
2242
2243 probdata = SCIPgetProbData(scip);
2244 assert(probdata != NULL);
2245 assert(probdata->vars != NULL);
2246
2247 nrsol = SCIPgetNSols(scip);
2248
2249 switch( SCIPgetStatus(scip) )
2250 {
2251 default:
2252 case SCIP_STATUS_UNKNOWN: /* the solving status is not yet known */
2253 gmoSolveStatSet(gmo, (int) gmoSolveStat_SystemErr);
2254 gmoModelStatSet(gmo, (int) gmoModelStat_ErrorNoSolution);
2255 nrsol = 0;
2256 break;
2257 case SCIP_STATUS_USERINTERRUPT: /* the user interrupted the solving process (by pressing Ctrl-C) */
2258 gmoSolveStatSet(gmo, (int) gmoSolveStat_User);
2259 gmoModelStatSet(gmo, nrsol > 0 ? (gmoNDisc(gmo) ? (int) gmoModelStat_Integer : (int) gmoModelStat_Feasible) : (int) gmoModelStat_NoSolutionReturned);
2260 break;
2261 case SCIP_STATUS_NODELIMIT: /* the solving process was interrupted because the node limit was reached */
2262 case SCIP_STATUS_STALLNODELIMIT: /* the solving process was interrupted because the node limit was reached */
2263 case SCIP_STATUS_TOTALNODELIMIT:
2264 gmoSolveStatSet(gmo, (int) gmoSolveStat_Iteration);
2265 gmoModelStatSet(gmo, nrsol > 0 ? (gmoNDisc(gmo) ? (int) gmoModelStat_Integer : (int) gmoModelStat_Feasible) : (int) gmoModelStat_NoSolutionReturned);
2266 break;
2267 case SCIP_STATUS_TIMELIMIT: /* the solving process was interrupted because the time limit was reached */
2268 case SCIP_STATUS_MEMLIMIT: /* the solving process was interrupted because the memory limit was reached */
2269 gmoSolveStatSet(gmo, (int) gmoSolveStat_Resource);
2270 gmoModelStatSet(gmo, nrsol > 0 ? (gmoNDisc(gmo) ? (int) gmoModelStat_Integer : (int) gmoModelStat_Feasible) : (int) gmoModelStat_NoSolutionReturned);
2271 break;
2272 case SCIP_STATUS_GAPLIMIT: /* the solving process was interrupted because the gap limit was reached */
2273 gmoSolveStatSet(gmo, (int) gmoSolveStat_Normal);
2274 gmoModelStatSet(gmo, nrsol > 0 ? (SCIPgetGap(scip) > 0.0 ? (gmoNDisc(gmo) ? (int) gmoModelStat_Integer : (int) gmoModelStat_Feasible) : (int) gmoModelStat_OptimalGlobal): (int) gmoModelStat_NoSolutionReturned);
2275 break;
2276 case SCIP_STATUS_SOLLIMIT: /* the solving process was interrupted because the solution limit was reached */
2277 case SCIP_STATUS_BESTSOLLIMIT: /* the solving process was interrupted because the solution improvement limit was reached */
2278 gmoSolveStatSet(gmo, (int) gmoSolveStat_Resource);
2279 gmoModelStatSet(gmo, nrsol > 0 ? (gmoNDisc(gmo) ? (int) gmoModelStat_Integer : (int) gmoModelStat_Feasible) : (int) gmoModelStat_NoSolutionReturned);
2280 break;
2281 case SCIP_STATUS_OPTIMAL: /* the problem was solved to optimality, an optimal solution is available */
2282 gmoSolveStatSet(gmo, (int) gmoSolveStat_Normal);
2283 gmoModelStatSet(gmo, (int) gmoModelStat_OptimalGlobal);
2284 break;
2285 case SCIP_STATUS_INFEASIBLE: /* the problem was proven to be infeasible */
2286 gmoSolveStatSet(gmo, (int) gmoSolveStat_Normal);
2287 gmoModelStatSet(gmo, (int) gmoModelStat_InfeasibleNoSolution);
2288 nrsol = 0;
2289 break;
2290 case SCIP_STATUS_UNBOUNDED: /* the problem was proven to be unbounded */
2291 gmoSolveStatSet(gmo, (int) gmoSolveStat_Normal);
2292 gmoModelStatSet(gmo, nrsol > 0 ? (int) gmoModelStat_Unbounded : (int) gmoModelStat_UnboundedNoSolution);
2293 break;
2294 case SCIP_STATUS_INFORUNBD: /* the problem was proven to be either infeasible or unbounded */
2295 gmoSolveStatSet(gmo, (int) gmoSolveStat_Normal);
2296 gmoModelStatSet(gmo, (int) gmoModelStat_NoSolutionReturned);
2297 nrsol = 0;
2298 break;
2299 }
2300
2301 if( SCIPgetStage(scip) == SCIP_STAGE_SOLVING || SCIPgetStage(scip) == SCIP_STAGE_SOLVED )
2302 dualbound = SCIPgetDualbound(scip);
2303 else
2304 dualbound = gmoValNA(gmo);
2305 gmoSetHeadnTail(gmo, (int) gmoTmipbest, dualbound);
2306 gmoSetHeadnTail(gmo, (int) gmoTmipnod, (double) SCIPgetNNodes(scip));
2307 gmoSetHeadnTail(gmo, (int) gmoHresused, SCIPgetSolvingTime(scip));
2308 gmoSetHeadnTail(gmo, (int) gmoHiterused, (double) SCIPgetNLPIterations(scip));
2309 gmoSetHeadnTail(gmo, (int) gmoHdomused, 0.0);
2310
2311 /* dump all solutions, if more than one found and parameter is set */
2312 if( nrsol > 1)
2313 {
2314 char* indexfilename;
2315
2316 SCIP_CALL( SCIPgetStringParam(scip, "gams/dumpsolutions", &indexfilename) );
2317 #ifndef WITH_GAMS
2318 if( indexfilename != NULL && indexfilename[0] )
2319 {
2320 char buffer[SCIP_MAXSTRLEN];
2321 gdxHandle_t gdx;
2322 int rc;
2323
2324 if( !gdxCreate(&gdx, buffer, sizeof(buffer)) )
2325 {
2326 SCIPerrorMessage("failed to load GDX I/O library: %s\n", buffer);
2327 return SCIP_OKAY;
2328 }
2329
2330 SCIPinfoMessage(scip, NULL, "\nDumping %d alternate solutions:\n", nrsol-1);
2331 /* create index GDX file */
2332 if( gdxOpenWrite(gdx, indexfilename, "SCIP DumpSolutions Index File", &rc) == 0 )
2333 {
2334 rc = gdxGetLastError(gdx);
2335 (void) gdxErrorStr(gdx, rc, buffer);
2336 SCIPerrorMessage("problem writing GDX file %s: %s\n", indexfilename, buffer);
2337 }
2338 else
2339 {
2340 gdxStrIndexPtrs_t keys;
2341 gdxStrIndex_t keysX;
2342 gdxValues_t vals;
2343 SCIP_Real* collev;
2344 int sloc;
2345 int i;
2346
2347 /* create index file */
2348 GDXSTRINDEXPTRS_INIT(keysX, keys);
2349 (void) gdxDataWriteStrStart(gdx, "index", "Dumpsolutions index", 1, (int) dt_set, 0);
2350 for( i = 1; i < nrsol; ++i)
2351 {
2352 (void) SCIPsnprintf(buffer, SCIP_MAXSTRLEN, "soln_scip_p%d.gdx", i);
2353 (void) gdxAddSetText(gdx, buffer, &sloc);
2354 (void) SCIPsnprintf(keys[0], GMS_SSSIZE, "file%d", i);
2355 vals[GMS_VAL_LEVEL] = sloc;
2356 (void) gdxDataWriteStr(gdx, (const char**)keys, vals);
2357 }
2358 (void) gdxDataWriteDone(gdx);
2359 (void) gdxClose(gdx);
2360
2361 SCIP_CALL( SCIPallocBufferArray(scip, &collev, gmoN(gmo)) );
2362
2363 /* create point files */
2364 for( i = 1; i < nrsol; ++i)
2365 {
2366 (void) SCIPsnprintf(buffer, SCIP_MAXSTRLEN, "soln_scip_p%d.gdx", i);
2367
2368 SCIP_CALL( SCIPgetSolVals(scip, SCIPgetSols(scip)[i], probdata->nvars, probdata->vars, collev) );
2369 (void) gmoSetVarL(gmo, collev);
2370 if( gmoUnloadSolutionGDX(gmo, buffer, 0, 1, 0) )
2371 {
2372 SCIPerrorMessage("Problems creating point file %s\n", buffer);
2373 }
2374 else
2375 {
2376 SCIPinfoMessage(scip, NULL, "Created point file %s\n", buffer);
2377 }
2378 }
2379
2380 SCIPfreeBufferArray(scip, &collev);
2381 }
2382
2383 (void) gdxFree(&gdx);
2384 }
2385
2386 SCIP_CALL( SCIPgetStringParam(scip, "gams/dumpsolutionsmerged", &indexfilename) );
2387 if( indexfilename != NULL && indexfilename[0] )
2388 {
2389 int solnvarsym;
2390
2391 if( gmoCheckSolPoolUEL(gmo, "soln_scip_p", &solnvarsym) )
2392 {
2393 SCIPerrorMessage("Solution pool scenario label 'soln_scip_p' contained in model dictionary. Cannot dump merged solutions pool.\n");
2394 }
2395 else
2396 {
2397 void* handle;
2398
2399 handle = gmoPrepareSolPoolMerge(gmo, indexfilename, nrsol-1, "soln_scip_p");
2400 if( handle != NULL )
2401 {
2402 SCIP_Real* collev;
2403 int k, i;
2404
2405 SCIP_CALL( SCIPallocBufferArray(scip, &collev, gmoN(gmo)) );
2406 for ( k = 0; k < solnvarsym; k++ )
2407 {
2408 gmoPrepareSolPoolNextSym(gmo, handle);
2409 for( i = 1; i < nrsol; ++i )
2410 {
2411 SCIP_CALL( SCIPgetSolVals(scip, SCIPgetSols(scip)[i], probdata->nvars, probdata->vars, collev) );
2412 (void) gmoSetVarL(gmo, collev);
2413 if( gmoUnloadSolPoolSolution (gmo, handle, i-1) )
2414 {
2415 SCIPerrorMessage("Problems unloading solution point %d symbol %d\n", i, k);
2416 }
2417 }
2418 }
2419 if( gmoFinalizeSolPoolMerge(gmo, handle) )
2420 {
2421 SCIPerrorMessage("Problems finalizing merged solution pool\n");
2422 }
2423 SCIPfreeBufferArray(scip, &collev);
2424 }
2425 else
2426 {
2427 SCIPerrorMessage("Problems preparing merged solution pool\n");
2428 }
2429 }
2430 }
2431 #endif
2432 }
2433
2434 /* pass best solution to GMO, if any */
2435 if( nrsol > 0 )
2436 {
2437 SCIP_SOL* sol;
2438 SCIP_Real* collev;
2439 SCIP_Real primalbound;
2440
2441 sol = SCIPgetBestSol(scip);
2442 assert(sol != NULL);
2443
2444 SCIP_CALL( SCIPallocBufferArray(scip, &collev, gmoN(gmo)) );
2445 SCIP_CALL( SCIPgetSolVals(scip, sol, probdata->nvars, probdata->vars, collev) );
2446
2447 #if GMOAPIVERSION < 12
2448 {
2449 SCIP_Real* lambda;
2450 int i;
2451
2452 SCIP_CALL( SCIPallocBufferArray(scip, &lambda, gmoM(gmo)) );
2453 for( i = 0; i < gmoM(gmo); ++i )
2454 lambda[i] = gmoValNA(gmo);
2455
2456 /* this also sets the gmoHobjval attribute to the level value of GAMS' objective variable */
2457 gmoSetSolution2(gmo, collev, lambda);
2458
2459 SCIPfreeBufferArray(scip, &lambda);
2460 }
2461 #else
2462 (void) gmoSetSolutionPrimal(gmo, collev);
2463 #endif
2464 primalbound = SCIPgetPrimalbound(scip);
2465
2466 SCIPfreeBufferArray(scip, &collev);
2467
2468 /* if we have an MINLP, check if best solution is really feasible and try to repair and check other solutions otherwise */
2469 if( gmoObjNLNZ(gmo) != 0 || gmoNLNZ(gmo) != 0 )
2470 {
2471 SCIP_Bool success;
2472
2473 /* check whether best solution is feasible in original problem */
2474 SCIP_CALL( SCIPcheckSolOrig(scip, sol, &success, FALSE, FALSE) );
2475
2476 /* look at all solutions, try to repair if not feasible, and keep a feasible one with best objective value */
2477 if( !success )
2478 {
2479 SCIP_Real mainsoltime;
2480 SCIP_CLOCK* resolveclock;
2481 int origmaxorigsol;
2482 int orignlpverblevel;
2483 int origmaxpresolrounds;
2484 char origsolvetracefile[SCIP_MAXSTRLEN];
2485 /* SCIP_Real origfeastol; */
2486 SCIP_Bool resolvenlp;
2487 SCIP_Real** solvals;
2488 SCIP_Real* objvals;
2489 int nsols;
2490 int s;
2491
2492 /* check whether we can run an NLP solver */
2493 SCIP_CALL( SCIPgetBoolParam(scip, "gams/resolvenlp", &resolvenlp) );
2494 if( resolvenlp && SCIPgetStage(scip) == SCIP_STAGE_SOLVING && !SCIPisNLPEnabled(scip) )
2495 {
2496 SCIPdebugMessage("NLP is disabled: cannot do resolves\n");
2497 resolvenlp = FALSE;
2498 }
2499 if( resolvenlp && SCIPgetNContVars(scip) == 0 )
2500 {
2501 SCIPdebugMessage("transformed SCIP problem has no continuous variables: cannot do resolves\n");
2502 resolvenlp = FALSE;
2503 }
2504 if( resolvenlp && SCIPgetNNlpis(scip) == 0 )
2505 {
2506 SCIPdebugMessage("no NLP solver: cannot do resolves\n");
2507 resolvenlp = FALSE;
2508 }
2509 if( resolvenlp && SCIPfindHeur(scip, "subnlp") == NULL )
2510 {
2511 SCIPdebugMessage("no NLP heuristic available: cannot do resolves\n");
2512 resolvenlp = FALSE;
2513 }
2514
2515 /* try SCIP solutions (limited by limits/maxsol or limits/maxorigsol) */
2516 nsols = SCIPgetNSols(scip);
2517 SCIP_CALL( SCIPallocBufferArray(scip, &solvals, nsols) );
2518 SCIP_CALL( SCIPallocBufferArray(scip, &objvals, nsols) );
2519 for( s = 0; s < nsols; ++s )
2520 {
2521 SCIP_CALL( SCIPallocBufferArray(scip, &solvals[s], gmoN(gmo)) ); /*lint !e866*/
2522 SCIP_CALL( SCIPgetSolVals(scip, SCIPgetSols(scip)[s], probdata->nvars, probdata->vars, solvals[s]) );
2523 objvals[s] = SCIPgetSolOrigObj(scip, SCIPgetSols(scip)[s]);
2524 }
2525
2526 /* adapt some parameter values for possible resolve's */
2527 if( resolvenlp )
2528 {
2529 char* tmp;
2530
2531 /* don't store solutions in original problem, so they don't get in a way when transforming for resolve */
2532 SCIP_CALL( SCIPgetIntParam(scip, "limits/maxorigsol", &origmaxorigsol) );
2533 SCIP_CALL( SCIPsetIntParam(scip, "limits/maxorigsol", 0) );
2534
2535 SCIP_CALL( SCIPgetIntParam(scip, "heuristics/subnlp/nlpverblevel", &orignlpverblevel) );
2536 SCIP_CALL( SCIPsetIntParam(scip, "heuristics/subnlp/nlpverblevel", 1) );
2537
2538 /* origfeastol = SCIPfeastol(scip); */
2539 /* SCIP_CALL( SCIPsetRealParam(scip, "numerics/feastol", origfeastol / 100.0) ); */
2540
2541 SCIP_CALL( SCIPgetIntParam(scip, "heuristics/subnlp/maxpresolverounds", &origmaxpresolrounds) );
2542 SCIP_CALL( SCIPsetIntParam(scip, "heuristics/subnlp/maxpresolverounds", 0) );
2543
2544 SCIP_CALL( SCIPgetStringParam(scip, "gams/solvetrace/file", &tmp) );
2545 strcpy(origsolvetracefile, tmp);
2546 SCIP_CALL( SCIPsetStringParam(scip, "gams/solvetrace/file", "") );
2547 }
2548
2549 SCIP_CALL( SCIPcreateClock(scip, &resolveclock) );
2550 SCIP_CALL( SCIPstartClock(scip, resolveclock) );
2551
2552 for( s = 0; s < nsols; ++s )
2553 {
2554 SCIPinfoMessage(scip, NULL, "Checking feasibility of solution #%.2d with reported objective value %.15e.\n", s, objvals[s]);
2555
2556 SCIP_CALL( checkAndRepairSol(scip, solvals[s], &objvals[s], resolvenlp, &success) );
2557
2558 if( success )
2559 break;
2560 }
2561
2562 SCIP_CALL( SCIPstopClock(scip, resolveclock) );
2563
2564 /* add time for checks and NLP resolves to reported solving time */
2565 mainsoltime = gmoGetHeadnTail(gmo, (int) gmoHresused);
2566 gmoSetHeadnTail(gmo, (int) gmoHresused, mainsoltime + SCIPgetClockTime(scip, resolveclock));
2567
2568 SCIP_CALL( SCIPfreeClock(scip, &resolveclock) );
2569
2570 if( success )
2571 {
2572 /* store updated solution in GMO */
2573 #if GMOAPIVERSION < 12
2574 {
2575 SCIP_Real* lambda;
2576 int i;
2577
2578 SCIP_CALL( SCIPallocBufferArray(scip, &lambda, gmoM(gmo)) );
2579 for( i = 0; i < gmoM(gmo); ++i )
2580 lambda[i] = gmoValNA(gmo);
2581
2582 /* this also sets the gmoHobjval attribute to the level value of GAMS' objective variable */
2583 gmoSetSolution2(gmo, solvals[s], lambda);
2584
2585 SCIPfreeBufferArray(scip, &lambda);
2586 }
2587 #else
2588 (void) gmoSetSolutionPrimal(gmo, solvals[s]);
2589 #endif
2590 /* update reevaluated objective value */
2591 objvals[s] = gmoGetHeadnTail(gmo, (int) gmoHobjval);
2592 primalbound = objvals[s];
2593
2594 SCIPinfoMessage(scip, NULL, "Solution #%.2d feasible. Reevaluated objective value = %.15e.\n", s, objvals[s]);
2595
2596 SCIPinfoMessage(scip, NULL, "\nStatus update:\n");
2597 SCIPinfoMessage(scip, NULL, "Solving Time (sec) : %.2f\n", gmoGetHeadnTail(gmo, (int) gmoHresused));
2598 SCIPinfoMessage(scip, NULL, "Primal Bound : %+.14e\n", objvals[s]);
2599 if( gmoGetHeadnTail(gmo, (int) gmoTmipbest) != gmoValNA(gmo) ) /*lint !e777*/
2600 {
2601 SCIPinfoMessage(scip, NULL, "Dual Bound : %+.14e\n", dualbound);
2602 SCIPinfoMessage(scip, NULL, "Gap : ");
2603
2604 if( SCIPisEQ(scip, objvals[s], dualbound) )
2605 SCIPinfoMessage(scip, NULL, "%.2f %%\n", 0.0);
2606 else if( SCIPisZero(scip, dualbound)
2607 || SCIPisZero(scip, objvals[s])
2608 || (dualbound == gmoValNA(gmo)) /*lint !e777*/
2609 || SCIPisInfinity(scip, REALABS(objvals[s]))
2610 || SCIPisInfinity(scip, REALABS(dualbound))
2611 || objvals[s] * dualbound < 0.0 )
2612 SCIPinfoMessage(scip, NULL, "infinite\n");
2613 else
2614 SCIPinfoMessage(scip, NULL, "%.2f %%\n", REALABS((objvals[s] - dualbound)/MIN(REALABS(dualbound),REALABS(objvals[s])))); /*lint !e666*/
2615 }
2616 }
2617 else
2618 {
2619 SCIPinfoMessage(scip, NULL, "None of %d SCIP solutions could be made feasible.\n", nsols);
2620 }
2621
2622 /* restore original parameter values */
2623 if( resolvenlp )
2624 {
2625 SCIP_CALL( SCIPsetIntParam(scip, "limits/maxorigsol", origmaxorigsol) ); /*lint !e644*/
2626 SCIP_CALL( SCIPsetIntParam(scip, "heuristics/subnlp/nlpverblevel", orignlpverblevel) ); /*lint !e644*/
2627 /* SCIP_CALL( SCIPsetRealParam(scip, "numerics/feastol", origfeastol) ); */
2628 SCIP_CALL( SCIPsetIntParam(scip, "heuristics/subnlp/maxpresolverounds", origmaxpresolrounds) ); /*lint !e644*/
2629 SCIP_CALL( SCIPsetStringParam(scip, "gams/solvetrace/file", origsolvetracefile) );
2630 }
2631
2632 for( s = 0; s < nsols; ++s )
2633 {
2634 SCIPfreeBufferArray(scip, &solvals[s]);
2635 }
2636 SCIPfreeBufferArray(scip, &solvals);
2637 SCIPfreeBufferArray(scip, &objvals);
2638 }
2639
2640 /* update model status */
2641 if( !success )
2642 {
2643 /* couldn't get a feasible solution, report intermediate infeasible */
2644 gmoModelStatSet(gmo, (int) gmoModelStat_InfeasibleIntermed);
2645 }
2646 else if( !SCIPisEQ(scip, primalbound, dualbound) )
2647 {
2648 /* feasible, but gap not closed, so only local optimum */
2649 gmoModelStatSet(gmo, gmoNDisc(gmo) ? (int) gmoModelStat_Integer : (int) gmoModelStat_Feasible);
2650 }
2651 else
2652 {
2653 /* feasible and gap closed, so solved globally */
2654 gmoModelStatSet(gmo, (int) gmoModelStat_OptimalGlobal);
2655 }
2656
2657 }
2658 }
2659
2660 if( gmoModelType(gmo) == (int) gmoProc_cns )
2661 switch( gmoModelStat(gmo) )
2662 {
2663 case gmoModelStat_OptimalGlobal:
2664 case gmoModelStat_OptimalLocal:
2665 case gmoModelStat_Feasible:
2666 case gmoModelStat_Integer:
2667 gmoModelStatSet(gmo, (int) gmoModelStat_Solved);
2668 gmoSetHeadnTail(gmo, gmoHobjval, 0.0);
2669 } /*lint !e744*/
2670
2671 return SCIP_OKAY;
2672 }
2673
2674 /*
2675 * Callback methods of reader
2676 */
2677
2678
2679 /** copy method for reader plugins (called when SCIP copies plugins) */
2680 #if 0
2681 static
2682 SCIP_DECL_READERCOPY(readerCopyGmo)
2683 { /*lint --e{715}*/
2684 SCIPerrorMessage("method of gmo reader not implemented yet\n");
2685 SCIPABORT(); /*lint --e{527}*/
2686
2687 return SCIP_OKAY;
2688 }
2689 #else
2690 #define readerCopyGmo NULL
2691 #endif
2692
2693 /** destructor of reader to free user data (called when SCIP is exiting) */
2694 #if 1
2695 static
SCIP_DECL_READERFREE(readerFreeGmo)2696 SCIP_DECL_READERFREE(readerFreeGmo)
2697 {
2698 SCIP_READERDATA* readerdata;
2699
2700 readerdata = SCIPreaderGetData(reader);
2701 assert(readerdata != NULL);
2702
2703 #if 0 /* TODO should do this if we created gmo */
2704 if( readerdata->gmo != NULL )
2705 {
2706 /* write solution file */
2707 gmoUnloadSolutionLegacy(readerdata->gmo);
2708
2709 gmoFree(&readerdata->gmo);
2710 gevFree(&readerdata->gev);
2711
2712 gmoLibraryUnload();
2713 gevLibraryUnload();
2714 }
2715 #endif
2716
2717 SCIPfreeMemory(scip, &readerdata);
2718
2719 return SCIP_OKAY;
2720 }
2721 #else
2722 #define readerFreeGmo NULL
2723 #endif
2724
2725
2726 /** problem reading method of reader */
2727 #if 1
2728 static
SCIP_DECL_READERREAD(readerReadGmo)2729 SCIP_DECL_READERREAD(readerReadGmo)
2730 {
2731 SCIP_READERDATA* readerdata;
2732 gmoHandle_t gmo;
2733 gevHandle_t gev;
2734 char buffer[1024];
2735
2736 *result = SCIP_DIDNOTRUN;
2737
2738 readerdata = SCIPreaderGetData(reader);
2739 assert(readerdata != NULL);
2740
2741 if( readerdata->gmo == NULL )
2742 {
2743 /* initialize GMO and GEV libraries */
2744 if( !gmoCreate(&readerdata->gmo, buffer, sizeof(buffer)) || !gevCreate(&readerdata->gev, buffer, sizeof(buffer)) )
2745 {
2746 SCIPerrorMessage(buffer);
2747 return SCIP_ERROR;
2748 }
2749
2750 gmo = readerdata->gmo;
2751 gev = readerdata->gev;
2752
2753 /* load control file */
2754 if( gevInitEnvironmentLegacy(gev, filename) )
2755 {
2756 SCIPerrorMessage("Could not load control file %s\n", filename);
2757 (void) gmoFree(&gmo);
2758 (void) gevFree(&gev);
2759 return SCIP_READERROR;
2760 }
2761
2762 if( gmoRegisterEnvironment(gmo, gev, buffer) )
2763 {
2764 SCIPerrorMessage("Error registering GAMS Environment: %s\n", buffer);
2765 (void) gmoFree(&gmo);
2766 (void) gevFree(&gev);
2767 return SCIP_ERROR;
2768 }
2769
2770 if( gmoLoadDataLegacy(gmo, buffer) )
2771 {
2772 SCIPerrorMessage("Could not load model data.\n");
2773 (void) gmoFree(&gmo);
2774 (void) gevFree(&gev);
2775 return SCIP_READERROR;
2776 }
2777 }
2778
2779 SCIP_CALL( SCIPcreateProblemReaderGmo(scip, readerdata->gmo, readerdata->indicatorfile, readerdata->mipstart) );
2780
2781 *result = SCIP_SUCCESS;
2782
2783 return SCIP_OKAY;
2784 }
2785 #else
2786 #define readerReadGmo NULL
2787 #endif
2788
2789
2790 #if 0
2791 /** problem writing method of reader */
2792 static
2793 SCIP_DECL_READERWRITE(readerWriteGmo)
2794 { /*lint --e{715}*/
2795 SCIPerrorMessage("method of gmo reader not implemented yet\n");
2796 SCIPABORT(); /*lint --e{527}*/
2797
2798 return SCIP_OKAY;
2799 }
2800 #else
2801 #define readerWriteGmo NULL
2802 #endif
2803
2804 /*
2805 * Constructs SCIP problem from the one in GMO.
2806 */
2807
2808 #define DIALOG_READGAMS_NAME "readgams"
2809 #define DIALOG_READGAMS_DESC "initializes SCIP problem to the one stored in a GAMS modeling object"
2810 #define DIALOG_READGAMS_ISSUBMENU FALSE
2811
2812 /** copy method for dialog plugins (called when SCIP copies plugins) */
2813 #if 0
2814 static
2815 SCIP_DECL_DIALOGCOPY(dialogCopyReadGams)
2816 { /*lint --e{715}*/
2817 SCIPerrorMessage("method of ReadGams dialog not implemented yet\n");
2818 SCIPABORT(); /*lint --e{527}*/
2819
2820 return SCIP_OKAY;
2821 }
2822 #else
2823 #define dialogCopyReadGams NULL
2824 #endif
2825
2826 /** destructor of dialog to free user data (called when the dialog is not captured anymore) */
2827 #if 0
2828 static
2829 SCIP_DECL_DIALOGFREE(dialogFreeReadGams)
2830 { /*lint --e{715}*/
2831 SCIPerrorMessage("method of ReadGams dialog not implemented yet\n");
2832 SCIPABORT(); /*lint --e{527}*/
2833
2834 return SCIP_OKAY;
2835 }
2836 #else
2837 #define dialogFreeReadGams NULL
2838 #endif
2839
2840 /** description output method of dialog */
2841 #if 0
2842 static
2843 SCIP_DECL_DIALOGDESC(dialogDescReadGams)
2844 { /*lint --e{715}*/
2845 SCIPerrorMessage("method of ReadGams dialog not implemented yet\n");
2846 SCIPABORT(); /*lint --e{527}*/
2847
2848 return SCIP_OKAY;
2849 }
2850 #else
2851 #define dialogDescReadGams NULL
2852 #endif
2853
2854
2855 /** execution method of dialog */
2856 static
SCIP_DECL_DIALOGEXEC(dialogExecReadGams)2857 SCIP_DECL_DIALOGEXEC(dialogExecReadGams)
2858 { /*lint --e{715}*/
2859 SCIP_READERDATA* readerdata;
2860
2861 assert(dialoghdlr != NULL);
2862 assert(dialog != NULL);
2863 assert(scip != NULL);
2864
2865 readerdata = (SCIP_READERDATA*)SCIPdialogGetData(dialog);
2866 assert(readerdata != NULL);
2867
2868 SCIP_CALL( SCIPdialoghdlrAddHistory(dialoghdlr, dialog, NULL, FALSE) );
2869
2870 if( readerdata->gmo == NULL )
2871 {
2872 SCIPerrorMessage("GMO not initialized, cannot setup GAMS problem\n");
2873 }
2874 else
2875 {
2876 /* free previous problem and solution data, if existing */
2877 SCIP_CALL( SCIPfreeProb(scip) );
2878
2879 SCIP_CALL( SCIPcreateProblemReaderGmo(scip, readerdata->gmo, readerdata->indicatorfile, readerdata->mipstart) );
2880 }
2881
2882 SCIPverbMessage(scip, SCIP_VERBLEVEL_NORMAL, NULL,
2883 "\noriginal problem has %d variables (%d bin, %d int, %d cont) and %d constraints\n",
2884 SCIPgetNOrigVars(scip), SCIPgetNOrigBinVars(scip), SCIPgetNOrigIntVars(scip), SCIPgetNOrigContVars(scip),
2885 SCIPgetNConss(scip));
2886
2887 *nextdialog = SCIPdialoghdlrGetRoot(dialoghdlr);
2888
2889 return SCIP_OKAY;
2890 }
2891
2892
2893 /*
2894 * Writing solution information to GMO dialog
2895 */
2896
2897 #define DIALOG_WRITEGAMSSOL_NAME "gamssol"
2898 #define DIALOG_WRITEGAMSSOL_DESC "writes solution information into GAMS Modeling Object"
2899 #define DIALOG_WRITEGAMSSOL_ISSUBMENU FALSE
2900
2901 /** copy method for dialog plugins (called when SCIP copies plugins) */
2902 #if 0
2903 static
2904 SCIP_DECL_DIALOGCOPY(dialogCopyWriteGamsSol)
2905 { /*lint --e{715}*/
2906 SCIPerrorMessage("method of WriteGamsSol dialog not implemented yet\n");
2907 SCIPABORT(); /*lint --e{527}*/
2908
2909 return SCIP_OKAY;
2910 }
2911 #else
2912 #define dialogCopyWriteGamsSol NULL
2913 #endif
2914
2915 /** destructor of dialog to free user data (called when the dialog is not captured anymore) */
2916 #if 0
2917 static
2918 SCIP_DECL_DIALOGFREE(dialogFreeWriteGamsSol)
2919 { /*lint --e{715}*/
2920 SCIPerrorMessage("method of WriteGamsSol dialog not implemented yet\n");
2921 SCIPABORT(); /*lint --e{527}*/
2922
2923 return SCIP_OKAY;
2924 }
2925 #else
2926 #define dialogFreeWriteGamsSol NULL
2927 #endif
2928
2929 /** description output method of dialog */
2930 #if 0
2931 static
2932 SCIP_DECL_DIALOGDESC(dialogDescWriteGamsSol)
2933 { /*lint --e{715}*/
2934 SCIPerrorMessage("method of WriteGamsSol dialog not implemented yet\n");
2935 SCIPABORT(); /*lint --e{527}*/
2936
2937 return SCIP_OKAY;
2938 }
2939 #else
2940 #define dialogDescWriteGamsSol NULL
2941 #endif
2942
2943
2944 /** execution method of dialog */
2945 static
SCIP_DECL_DIALOGEXEC(dialogExecWriteGamsSol)2946 SCIP_DECL_DIALOGEXEC(dialogExecWriteGamsSol)
2947 { /*lint --e{715}*/
2948 SCIP_READERDATA* readerdata;
2949
2950 SCIP_CALL( SCIPdialoghdlrAddHistory(dialoghdlr, dialog, NULL, FALSE) );
2951
2952 readerdata = (SCIP_READERDATA*) SCIPdialogGetData(dialog);
2953 assert(readerdata != NULL);
2954
2955 SCIP_CALL( writeGmoSolution(scip, readerdata->gmo) );
2956
2957 *nextdialog = SCIPdialoghdlrGetRoot(dialoghdlr);
2958
2959 return SCIP_OKAY;
2960 }
2961
2962
2963 /*
2964 * Loading GAMS and user option file dialog
2965 */
2966
2967 #define DIALOG_SETTINGSLOADGAMS_NAME "loadgams"
2968 #define DIALOG_SETTINGSLOADGAMS_DESC "loads GAMS settings and SCIP option file specified in GAMS model"
2969 #define DIALOG_SETTINGSLOADGAMS_ISSUBMENU FALSE
2970
2971 /** copy method for dialog plugins (called when SCIP copies plugins) */
2972 #if 0
2973 static
2974 SCIP_DECL_DIALOGCOPY(dialogCopySettingsLoadGams)
2975 { /*lint --e{715}*/
2976 SCIPerrorMessage("method of SettingsLoadGams dialog not implemented yet\n");
2977 SCIPABORT(); /*lint --e{527}*/
2978
2979 return SCIP_OKAY;
2980 }
2981 #else
2982 #define dialogCopySettingsLoadGams NULL
2983 #endif
2984
2985 /** destructor of dialog to free user data (called when the dialog is not captured anymore) */
2986 #if 0
2987 static
2988 SCIP_DECL_DIALOGFREE(dialogFreeSettingsLoadGams)
2989 { /*lint --e{715}*/
2990 SCIPerrorMessage("method of SettingsLoadGams dialog not implemented yet\n");
2991 SCIPABORT(); /*lint --e{527}*/
2992
2993 return SCIP_OKAY;
2994 }
2995 #else
2996 #define dialogFreeSettingsLoadGams NULL
2997 #endif
2998
2999 /** description output method of dialog */
3000 #if 0
3001 static
3002 SCIP_DECL_DIALOGDESC(dialogDescSettingsLoadGams)
3003 { /*lint --e{715}*/
3004 SCIPerrorMessage("method of SettingsLoadGams dialog not implemented yet\n");
3005 SCIPABORT(); /*lint --e{527}*/
3006
3007 return SCIP_OKAY;
3008 }
3009 #else
3010 #define dialogDescSettingsLoadGams NULL
3011 #endif
3012
3013
3014 /** execution method of dialog */
3015 static
SCIP_DECL_DIALOGEXEC(dialogExecSettingsLoadGams)3016 SCIP_DECL_DIALOGEXEC(dialogExecSettingsLoadGams)
3017 { /*lint --e{715}*/
3018 assert(scip != NULL);
3019 assert(dialog != NULL);
3020 assert(dialoghdlr != NULL);
3021
3022 SCIP_CALL( SCIPdialoghdlrAddHistory(dialoghdlr, dialog, NULL, FALSE) );
3023
3024 SCIP_CALL( SCIPreadParamsReaderGmo(scip) );
3025
3026 *nextdialog = SCIPdialoghdlrGetRoot(dialoghdlr);
3027
3028 return SCIP_OKAY;
3029 }
3030
3031 /*
3032 * reader specific interface methods
3033 */
3034
3035 /** includes the gmo file reader in SCIP */
SCIPincludeReaderGmo(SCIP * scip)3036 SCIP_RETCODE SCIPincludeReaderGmo(
3037 SCIP* scip /**< SCIP data structure */
3038 )
3039 {
3040 SCIP_READERDATA* readerdata;
3041 SCIP_DIALOG* dialog;
3042 SCIP_DIALOG* parentdialog;
3043
3044 /* create gmo reader data */
3045 SCIP_CALL( SCIPallocMemory(scip, &readerdata) );
3046 BMSclearMemory(readerdata);
3047
3048 /* include gmo reader */
3049 SCIP_CALL( SCIPincludeReader(scip, READER_NAME, READER_DESC, READER_EXTENSION,
3050 readerCopyGmo,
3051 readerFreeGmo, readerReadGmo, readerWriteGmo,
3052 readerdata) );
3053
3054 SCIP_CALL( SCIPaddStringParam(scip, "gams/dumpsolutions",
3055 "name of solutions index gdx file for writing all alternate solutions",
3056 NULL, FALSE, "", NULL, NULL) );
3057
3058 SCIP_CALL( SCIPaddStringParam(scip, "gams/dumpsolutionsmerged",
3059 "name of gdx file for writing all alternate solutions into a single file",
3060 NULL, FALSE, "", NULL, NULL) );
3061
3062 SCIP_CALL( SCIPaddBoolParam(scip, "gams/resolvenlp",
3063 "whether to resolve MINLP with fixed discrete variables if best solution violates some constraints",
3064 NULL, FALSE, TRUE, NULL, NULL) );
3065
3066 SCIP_CALL( SCIPaddIntParam(scip, "gams/mipstart",
3067 "how to handle initial variable levels",
3068 &readerdata->mipstart, FALSE, 2, 0, 4, NULL, NULL) );
3069
3070 SCIP_CALL( SCIPaddStringParam(scip, "gams/indicatorfile",
3071 "name of GAMS options file that contains definitions on indicators",
3072 &readerdata->indicatorfile, FALSE, "", NULL, NULL) );
3073
3074 /* get parent dialog "write" */
3075 if( SCIPdialogFindEntry(SCIPgetRootDialog(scip), "write", &parentdialog) != 1 )
3076 {
3077 SCIPerrorMessage("sub menu \"write\" not found\n");
3078 return SCIP_PLUGINNOTFOUND;
3079 }
3080 assert(parentdialog != NULL);
3081
3082 /* create, include, and release dialog */
3083 if( !SCIPdialogHasEntry(SCIPgetRootDialog(scip), DIALOG_READGAMS_NAME) )
3084 {
3085 SCIP_CALL( SCIPincludeDialog(scip, &dialog,
3086 dialogCopyReadGams, dialogExecReadGams, dialogDescReadGams, dialogFreeReadGams,
3087 DIALOG_READGAMS_NAME, DIALOG_READGAMS_DESC, DIALOG_READGAMS_ISSUBMENU, (SCIP_DIALOGDATA*)readerdata) );
3088 SCIP_CALL( SCIPaddDialogEntry(scip, SCIPgetRootDialog(scip), dialog) );
3089 SCIP_CALL( SCIPreleaseDialog(scip, &dialog) );
3090 }
3091
3092
3093 /* create, include, and release dialog */
3094 if( !SCIPdialogHasEntry(parentdialog, DIALOG_WRITEGAMSSOL_NAME) )
3095 {
3096 SCIP_CALL( SCIPincludeDialog(scip, &dialog,
3097 dialogCopyWriteGamsSol, dialogExecWriteGamsSol, dialogDescWriteGamsSol, dialogFreeWriteGamsSol,
3098 DIALOG_WRITEGAMSSOL_NAME, DIALOG_WRITEGAMSSOL_DESC, DIALOG_WRITEGAMSSOL_ISSUBMENU, (SCIP_DIALOGDATA*)readerdata) );
3099 SCIP_CALL( SCIPaddDialogEntry(scip, parentdialog, dialog) );
3100 SCIP_CALL( SCIPreleaseDialog(scip, &dialog) );
3101 }
3102
3103
3104 /* get parent dialog "set" */
3105 if( SCIPdialogFindEntry(SCIPgetRootDialog(scip), "set", &parentdialog) != 1 )
3106 {
3107 SCIPerrorMessage("sub menu \"set\" not found\n");
3108 return SCIP_PLUGINNOTFOUND;
3109 }
3110 assert(parentdialog != NULL);
3111
3112 /* create, include, and release dialog */
3113 if( !SCIPdialogHasEntry(parentdialog, DIALOG_SETTINGSLOADGAMS_NAME) )
3114 {
3115 SCIP_CALL( SCIPincludeDialog(scip, &dialog,
3116 dialogCopySettingsLoadGams, dialogExecSettingsLoadGams, dialogDescSettingsLoadGams, dialogFreeSettingsLoadGams,
3117 DIALOG_SETTINGSLOADGAMS_NAME, DIALOG_SETTINGSLOADGAMS_DESC, DIALOG_SETTINGSLOADGAMS_ISSUBMENU, NULL) );
3118 SCIP_CALL( SCIPaddDialogEntry(scip, parentdialog, dialog) );
3119 SCIP_CALL( SCIPreleaseDialog(scip, &dialog) );
3120 }
3121
3122 /* set gams */
3123 if( !SCIPdialogHasEntry(parentdialog, "gams") )
3124 {
3125 SCIP_CALL( SCIPincludeDialog(scip, &dialog,
3126 NULL, SCIPdialogExecMenu, NULL, NULL,
3127 "gams", "change parameters for GAMS interface", TRUE, NULL) );
3128 SCIP_CALL( SCIPaddDialogEntry(scip, parentdialog, dialog) );
3129 SCIP_CALL( SCIPreleaseDialog(scip, &dialog) );
3130 }
3131
3132 return SCIP_OKAY;
3133 }
3134
3135 /** sets the GMO object to use in reader
3136 * If GMO is set in reader, then reader does not read from file when executed, but sets up problem from GMO
3137 */
SCIPsetGMOReaderGmo(SCIP * scip,gmoRec_t * gmo)3138 void SCIPsetGMOReaderGmo(
3139 SCIP* scip, /**< SCIP data structure */
3140 gmoRec_t* gmo /**< GMO object, or NULL to reset to default behaviour */
3141 )
3142 {
3143 SCIP_READER* reader;
3144 SCIP_READERDATA* readerdata;
3145
3146 assert(scip != NULL);
3147 assert(SCIPgetStage(scip) == SCIP_STAGE_INIT);
3148
3149 reader = SCIPfindReader(scip, READER_NAME);
3150 assert(reader != NULL);
3151
3152 readerdata = SCIPreaderGetData(reader);
3153 assert(readerdata != NULL);
3154
3155 readerdata->gmo = gmo;
3156 readerdata->gev = gmo != NULL ? (gevHandle_t)gmoEnvironment(readerdata->gmo) : NULL;
3157 }
3158
3159 /** passes GAMS options to SCIP and initiates reading of user options file, if given in GMO */
SCIPreadParamsReaderGmo(SCIP * scip)3160 SCIP_RETCODE SCIPreadParamsReaderGmo(
3161 SCIP* scip /**< SCIP data structure */
3162 )
3163 {
3164 SCIP_READER* reader;
3165 SCIP_READERDATA* readerdata;
3166 gmoHandle_t gmo;
3167 gevHandle_t gev;
3168
3169 assert(scip != NULL);
3170
3171 reader = SCIPfindReader(scip, READER_NAME);
3172 assert(reader != NULL);
3173
3174 readerdata = SCIPreaderGetData(reader);
3175 assert(readerdata != NULL);
3176 assert(readerdata->gmo != NULL);
3177 assert(readerdata->gev != NULL);
3178
3179 gmo = readerdata->gmo;
3180 gev = readerdata->gev;
3181
3182 if( gevGetIntOpt(gev, gevNodeLim) > 0 )
3183 {
3184 SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", (long long)gevGetIntOpt(gev, gevNodeLim)) );
3185 }
3186 SCIP_CALL( SCIPsetRealParam(scip, "limits/time", gevGetDblOpt(gev, gevResLim)) );
3187 SCIP_CALL( SCIPsetRealParam(scip, "limits/gap", gevGetDblOpt(gev, gevOptCR)) );
3188 if( !SCIPisInfinity(scip, gevGetDblOpt(gev, gevOptCA)) )
3189 {
3190 SCIP_CALL( SCIPsetRealParam(scip, "limits/absgap", gevGetDblOpt(gev, gevOptCA)) );
3191 }
3192 else
3193 {
3194 SCIPwarningMessage(scip, "Value for optca = %g >= value for infinity. Setting solution limit to 1 instead.\n", gevGetDblOpt(gev, gevOptCA));
3195 SCIP_CALL( SCIPsetIntParam(scip, "limits/solutions", 1) );
3196 }
3197 if( gevGetDblOpt(gev, gevWorkSpace) > 0.0 )
3198 {
3199 SCIP_CALL( SCIPsetRealParam(scip, "limits/memory", gevGetDblOpt(gev, gevWorkSpace)) );
3200 }
3201 SCIP_CALL( SCIPsetIntParam(scip, "lp/threads", gevThreads(gev)) );
3202
3203 /* if log is not kept, then can also set SCIP verblevel to 0 */
3204 if( gevGetIntOpt(gev, gevLogOption) == 0 )
3205 {
3206 SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 0) );
3207 }
3208
3209 #ifdef _WIN32
3210 if( !gevGetIntOpt(gev, gevIDEFlag) )
3211 {
3212 SCIP_CALL( SCIPsetIntParam(scip, "display/width", 80) );
3213 SCIP_CALL( SCIPsetIntParam(scip, "display/lpavgiterations/active", 0) );
3214 SCIP_CALL( SCIPsetIntParam(scip, "display/maxdepth/active", 0) );
3215 SCIP_CALL( SCIPsetIntParam(scip, "display/time/active", 2) );
3216 }
3217 #endif
3218
3219 /* enable column on number of branching on nonlinear variables, if any */
3220 if( gmoNLNZ(gmo) > 0 || (gmoObjStyle(gmo) == (int) gmoObjType_Fun && gmoObjNLNZ(gmo) > 0) )
3221 {
3222 /* enable column on number of branching on continuous variables */
3223 SCIP_CALL( SCIPsetIntParam(scip, "display/nexternbranchcands/active", 2) );
3224 }
3225 /* make sure column on number of branching on fractional variables is shown, if any */
3226 if( gmoNDisc(gmo) > 0 )
3227 {
3228 SCIP_CALL( SCIPsetIntParam(scip, "display/nfrac/active", 2) );
3229 }
3230
3231 /* don't print reason why start solution is infeasible, per default */
3232 SCIP_CALL( SCIPsetBoolParam(scip, "misc/printreason", FALSE) );
3233
3234 if( gmoOptFile(gmo) > 0 )
3235 {
3236 char optfilename[1024];
3237 SCIP_RETCODE ret;
3238
3239 (void) gmoNameOptFile(gmo, optfilename);
3240 SCIPinfoMessage(scip, NULL, "\nreading option file %s\n", optfilename);
3241 ret = SCIPreadParams(scip, optfilename);
3242 if( ret != SCIP_OKAY )
3243 {
3244 SCIPwarningMessage(scip, "Reading of optionfile %s failed with SCIP return code <%d>!\n", optfilename, ret);
3245 }
3246 }
3247
3248 return SCIP_OKAY;
3249 }
3250