1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file reader_nl.cpp
17 * @ingroup DEFPLUGINS_READER
18 * @brief AMPL .nl file reader
19 * @author Stefan Vigerske
20 *
21 * For documentation on ampl::mp, see https://ampl.github.io and https://www.zverovich.net/2014/09/19/reading-nl-files.html.
22 * For documentation on .nl files, see https://ampl.com/REFS/hooking2.pdf.
23 */
24
25 /*--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
26
27 #include <string>
28 #include <vector>
29 #include <map>
30
31 #include "scip/reader_nl.h"
32 #include "scip/cons_linear.h"
33 #include "scip/cons_nonlinear.h"
34 #include "scip/cons_sos1.h"
35 #include "scip/cons_sos2.h"
36 #include "scip/expr_var.h"
37 #include "scip/expr_value.h"
38 #include "scip/expr_sum.h"
39 #include "scip/expr_product.h"
40 #include "scip/expr_pow.h"
41 #include "scip/expr_log.h"
42 #include "scip/expr_exp.h"
43 #include "scip/expr_trig.h"
44 #include "scip/expr_abs.h"
45
46 // disable -Wshadow warnings for upcoming includes of AMPL/MP
47 // disable -Wimplicit-fallthrough as I don't want to maintain extra comments in AMPL/MP code to suppress these
48 #ifdef __GNUC__
49 #pragma GCC diagnostic ignored "-Wshadow"
50 #if __GNUC__ >= 7
51 #pragma GCC diagnostic ignored "-Wimplicit-fallthrough"
52 #endif
53 #endif
54
55 #include "mp/nl-reader.h"
56
57 #define READER_NAME "nlreader"
58 #define READER_DESC "AMPL .nl file reader"
59 #define READER_EXTENSION "nl"
60
61 // a variant of SCIP_CALL that throws a std::logic_error if not SCIP_OKAY
62 // (using cast to long long to work around issues with old MSVC)
63 #define SCIP_CALL_THROW(x) \
64 do \
65 { \
66 SCIP_RETCODE throw_retcode; \
67 if( ((throw_retcode) = (x)) != SCIP_OKAY ) \
68 throw std::logic_error("Error <" + std::to_string((long long)throw_retcode) + "> in function call"); \
69 } \
70 while( false )
71
72 /*
73 * Data structures
74 */
75
76 /// problem data stored in SCIP
77 struct SCIP_ProbData
78 {
79 char* filenamestub; /**< name of input file, without .nl extension; array is long enough to hold 5 extra chars */
80 int filenamestublen; /**< length of filenamestub string */
81
82 int amplopts[mp::MAX_AMPL_OPTIONS]; /**< AMPL options from .nl header */
83 int namplopts; /**< number of AMPL options from .nl header */
84
85 SCIP_VAR** vars; /**< variables in the order given by AMPL */
86 int nvars; /**< number of variables */
87
88 SCIP_CONS** conss; /**< constraints in the order given by AMPL */
89 int nconss; /**< number of constraints */
90
91 SCIP_Bool islp; /**< whether problem is an LP (only linear constraints, only continuous vars) */
92 };
93
94 /*
95 * Local methods
96 */
97
98 // forward declaration
99 static SCIP_DECL_PROBDELORIG(probdataDelOrigNl);
100
101 /// implementation of AMPL/MPs NLHandler that constructs a SCIP problem while a .nl file is read
102 class AMPLProblemHandler : public mp::NLHandler<AMPLProblemHandler, SCIP_EXPR*>
103 {
104 private:
105 SCIP* scip;
106 SCIP_PROBDATA* probdata;
107
108 // variable expressions corresponding to nonlinear variables
109 // created in OnHeader() and released in destructor
110 // for reuse of var-expressions in OnVariableRef()
111 std::vector<SCIP_EXPR*> varexprs;
112
113 // linear parts for nonlinear constraints
114 // first collect and then add to constraints in EndInput()
115 std::vector<std::vector<std::pair<SCIP_Real, SCIP_VAR*> > > nlconslin;
116
117 // expression that represents a nonlinear objective function
118 // used to create a corresponding constraint in EndInput(), unless NULL
119 SCIP_EXPR* objexpr;
120
121 // common expressions (defined variables from statements like "var xsqr = x^2;" in an AMPL model)
122 // they are constructed by BeginCommonExpr/EndCommonExpr below and are referenced by index in OnCommonExprRef
123 std::vector<SCIP_EXPR*> commonexprs;
124
125 // collect expressions that need to be released eventually
126 // this are all expression that are returned to the AMPL/MP code in AMPLProblemHandler::OnXyz() functions
127 // they need to be released exactly once, but after they are used in another expression or a constraint
128 // as AMPL/MP may reuse expressions (common subexpressions), we don't release an expression when it is used
129 // as a child or when constructing a constraint, but first collect them all and then release in destructor
130 // alternatively, one could encapsulate SCIP_EXPR* into a small class that handles proper reference counting
131 std::vector<SCIP_EXPR*> exprstorelease;
132
133 // SOS constraints
134 // collected while handling suffixes in SuffixHandler
135 // sosvars maps the SOS index (can be negative) to the indices of the variables in the SOS
136 // sosweights gives for each variable its weight in the SOS it appears in (if any)
137 std::map<int, std::vector<int> > sosvars;
138 std::vector<int> sosweights;
139
140 // initial solution, if any
141 SCIP_SOL* initsol;
142
143 // opened files with column/variable and row/constraint names, or NULL
144 fmt::File* colfile;
145 fmt::File* rowfile;
146
147 // get name from names strings, if possible
148 // returns whether a name has been stored
149 bool nextName(
150 const char*& namesbegin, /**< current pointer into names string, or NULL */
151 const char* namesend, /**< pointer to end of names string */
152 char* name /**< buffer to store name, should have length SCIP_MAXSTRLEN */
153 )
154 {
155 if( namesbegin == NULL )
156 return false;
157
158 // copy namesbegin into name until newline or namesend
159 // updates namesbegin
160 int nchars = 0;
161 while( namesbegin != namesend )
162 {
163 if( nchars == SCIP_MAXSTRLEN )
164 {
165 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "name too long when parsing names file");
166 // do no longer read names from this string (something seems awkward)
167 namesbegin = NULL;
168 return false;
169 }
170 if( *namesbegin == '\n' )
171 {
172 *name = '\0';
173 ++namesbegin;
174 return true;
175 }
176 *(name++) = *(namesbegin++);
177 ++nchars;
178 }
179
180 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "missing newline when parsing names file");
181 return false;
182 }
183
184 public:
185 /// constructor
186 ///
187 /// initializes SCIP problem and problem data
188 AMPLProblemHandler(
189 SCIP* scip_, ///< SCIP data structure
190 const char* filename ///< name of .nl file that is read
191 )
192 : scip(scip_),
193 probdata(NULL),
194 objexpr(NULL),
195 initsol(NULL),
196 colfile(NULL),
197 rowfile(NULL)
198 {
199 assert(scip != NULL);
200 assert(filename != NULL);
201
202 SCIP_CALL_THROW( SCIPallocClearMemory(scip, &probdata) );
203
204 /* get name of input file without file extension (if any) */
205 const char* extstart = strrchr(const_cast<char*>(filename), '.');
206 if( extstart != NULL )
207 probdata->filenamestublen = extstart - filename;
208 else
209 probdata->filenamestublen = strlen(filename);
210 assert(probdata->filenamestublen > 0);
211 SCIP_CALL_THROW( SCIPallocBlockMemoryArray(scip, &probdata->filenamestub, probdata->filenamestublen + 5) );
212 memcpy(probdata->filenamestub, filename, probdata->filenamestublen);
213 probdata->filenamestub[probdata->filenamestublen] = '\0';
214
215 /* derive probname from name of input file without path and extension */
216 const char* probname = strrchr(probdata->filenamestub, '/');
217 if( probname == NULL )
218 probname = probdata->filenamestub;
219 else
220 ++probname;
221
222 // initialize empty SCIP problem
223 SCIP_CALL_THROW( SCIPcreateProb(scip, probname, probdataDelOrigNl, NULL, NULL, NULL, NULL, NULL, probdata) );
224
225 // try to open files with variable and constraint names
226 // temporarily add ".col" and ".row", respectively, to filenamestub
227 try
228 {
229 probdata->filenamestub[probdata->filenamestublen] = '.';
230 probdata->filenamestub[probdata->filenamestublen+1] = 'c';
231 probdata->filenamestub[probdata->filenamestublen+2] = 'o';
232 probdata->filenamestub[probdata->filenamestublen+3] = 'l';
233 probdata->filenamestub[probdata->filenamestublen+4] = '\0';
234 colfile = new fmt::File(probdata->filenamestub, fmt::File::RDONLY);
235
236 probdata->filenamestub[probdata->filenamestublen+1] = 'r';
237 probdata->filenamestub[probdata->filenamestublen+3] = 'w';
238 rowfile = new fmt::File(probdata->filenamestub, fmt::File::RDONLY);
239 }
240 catch( const fmt::SystemError& e )
241 {
242 // probably a file open error, probably because file not found
243 // ignore, we can make up our own names
244 }
245 probdata->filenamestub[probdata->filenamestublen] = '\0';
246 }
247
248 /// destructor
249 ///
250 /// only asserts that cleanup() has been called, as we cannot throw an exception or return a SCIP_RETCODE here
251 ~AMPLProblemHandler()
252 {
253 // exprs and linear constraint arrays should have been cleared up in cleanup()
254 assert(varexprs.empty());
255 assert(exprstorelease.empty());
256
257 delete colfile;
258 delete rowfile;
259 }
260
261 /// process header of .nl files
262 ///
263 /// create and add variables, allocate constraints
264 void OnHeader(
265 const mp::NLHeader& h ///< header data
266 )
267 {
268 char name[SCIP_MAXSTRLEN];
269 int nnlvars;
270
271 assert(probdata->vars == NULL);
272 assert(probdata->conss == NULL);
273
274 probdata->namplopts = h.num_ampl_options;
275 BMScopyMemoryArray(probdata->amplopts, h.ampl_options, h.num_ampl_options);
276
277 // read variable and constraint names from file, if available, into memory
278 // if not available, we will get varnamesbegin==NULL and consnamesbegin==NULL
279 mp::MemoryMappedFile<> mapped_colfile;
280 if( colfile != NULL )
281 mapped_colfile.map(*colfile, "colfile");
282 const char* varnamesbegin = mapped_colfile.start();
283 const char* varnamesend = mapped_colfile.start() + mapped_colfile.size();
284
285 mp::MemoryMappedFile<> mapped_rowfile;
286 if( rowfile != NULL )
287 mapped_rowfile.map(*rowfile, "rowfile");
288 const char* consnamesbegin = mapped_rowfile.start();
289 const char* consnamesend = mapped_rowfile.start() + mapped_rowfile.size();
290
291 probdata->nvars = h.num_vars;
292 SCIP_CALL_THROW( SCIPallocBlockMemoryArray(scip, &probdata->vars, probdata->nvars) );
293
294 // number of nonlinear variables
295 nnlvars = MAX(h.num_nl_vars_in_cons, h.num_nl_vars_in_objs);
296 varexprs.resize(nnlvars);
297
298 // create variables
299 // create variable expressions for nonlinear variables
300 for( int i = 0; i < h.num_vars; ++i )
301 {
302 SCIP_VARTYPE vartype;
303 // Nonlinear variables in both constraints and objective
304 if( i < h.num_nl_vars_in_both - h.num_nl_integer_vars_in_both )
305 vartype = SCIP_VARTYPE_CONTINUOUS;
306 else if( i < h.num_nl_vars_in_both )
307 vartype = SCIP_VARTYPE_INTEGER;
308 // Nonlinear variables in constraints
309 else if( i < h.num_nl_vars_in_cons - h.num_nl_integer_vars_in_cons )
310 vartype = SCIP_VARTYPE_CONTINUOUS;
311 else if( i < h.num_nl_vars_in_cons )
312 vartype = SCIP_VARTYPE_INTEGER;
313 // Nonlinear variables in objective
314 else if( i < h.num_nl_vars_in_objs - h.num_nl_integer_vars_in_objs )
315 vartype = SCIP_VARTYPE_CONTINUOUS;
316 else if( i < h.num_nl_vars_in_objs )
317 vartype = SCIP_VARTYPE_INTEGER;
318 // Linear variables
319 else if( i < h.num_vars - h.num_linear_binary_vars - h.num_linear_integer_vars )
320 vartype = SCIP_VARTYPE_CONTINUOUS;
321 else if( i < h.num_vars - h.num_linear_integer_vars )
322 vartype = SCIP_VARTYPE_BINARY;
323 else
324 vartype = SCIP_VARTYPE_INTEGER;
325
326 if( !nextName(varnamesbegin, varnamesend, name) )
327 {
328 // make up name if no names file or could not be read
329 switch( vartype )
330 {
331 case SCIP_VARTYPE_BINARY :
332 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "b%d", i);
333 break;
334 case SCIP_VARTYPE_INTEGER :
335 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "i%d", i);
336 break;
337 case SCIP_VARTYPE_CONTINUOUS :
338 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x%d", i);
339 break;
340 default:
341 SCIPABORT();
342 break;
343 }
344 }
345
346 SCIP_CALL_THROW( SCIPcreateVarBasic(scip, &probdata->vars[i], name,
347 vartype == SCIP_VARTYPE_BINARY ? 0.0 : -SCIPinfinity(scip),
348 vartype == SCIP_VARTYPE_BINARY ? 1.0 : SCIPinfinity(scip),
349 0.0, vartype) );
350 SCIP_CALL_THROW( SCIPaddVar(scip, probdata->vars[i]) );
351
352 if( i < nnlvars )
353 {
354 SCIP_CALL_THROW( SCIPcreateExprVar(scip, &varexprs[i], probdata->vars[i], NULL, NULL) );
355 }
356 }
357
358 // alloc some space for constraints
359 probdata->nconss = h.num_algebraic_cons;
360 SCIP_CALL_THROW( SCIPallocBlockMemoryArray(scip, &probdata->conss, probdata->nconss) );
361 nlconslin.resize(h.num_nl_cons);
362
363 // create empty nonlinear constraints
364 // use expression == 0, because nonlinear constraint don't like to be without an expression
365 SCIP_EXPR* dummyexpr;
366 SCIP_CALL_THROW( SCIPcreateExprValue(scip, &dummyexpr, 0.0, NULL, NULL) );
367 for( int i = 0; i < h.num_nl_cons; ++i )
368 {
369 // make up name if no names file or could not be read
370 if( !nextName(consnamesbegin, consnamesend, name) )
371 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlc%d", i);
372
373 SCIP_CALL_THROW( SCIPcreateConsBasicNonlinear(scip, &probdata->conss[i], name, dummyexpr, -SCIPinfinity(scip), SCIPinfinity(scip)) );
374 }
375 SCIP_CALL_THROW( SCIPreleaseExpr(scip, &dummyexpr) );
376
377 // create empty linear constraints
378 for( int i = h.num_nl_cons; i < h.num_algebraic_cons; ++i )
379 {
380 if( !nextName(consnamesbegin, consnamesend, name) )
381 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "lc%d", i);
382 SCIP_CALL_THROW( SCIPcreateConsBasicLinear(scip, &probdata->conss[i], name, 0, NULL, NULL, -SCIPinfinity(scip), SCIPinfinity(scip)) );
383 }
384
385 if( h.num_nl_cons == 0 && h.num_integer_vars() == 0 )
386 probdata->islp = true;
387
388 // alloc space for common expressions
389 commonexprs.resize(h.num_common_exprs());
390 }
391
392 /// receive notification of a number in a nonlinear expression
393 SCIP_EXPR* OnNumber(
394 double value ///< value
395 )
396 {
397 SCIP_EXPR* expr;
398
399 SCIP_CALL_THROW( SCIPcreateExprValue(scip, &expr, value, NULL, NULL) );
400
401 // remember that we have to release this expr
402 exprstorelease.push_back(expr);
403
404 return expr;
405 }
406
407 /// receive notification of a variable reference in a nonlinear expression
408 SCIP_EXPR* OnVariableRef(
409 int variableIndex ///< AMPL index of variable
410 )
411 {
412 assert(variableIndex >= 0);
413 assert(variableIndex < (int)varexprs.size());
414 assert(varexprs[variableIndex] != NULL);
415
416 return varexprs[variableIndex];
417 }
418
419 /// receive notification of a unary expression
420 SCIP_EXPR* OnUnary(
421 mp::expr::Kind kind, ///< expression operator
422 SCIP_EXPR* child ///< argument
423 )
424 {
425 SCIP_EXPR* expr;
426
427 assert(child != NULL);
428
429 switch( kind )
430 {
431 case mp::expr::MINUS:
432 {
433 SCIP_Real minusone = -1.0;
434 SCIP_CALL_THROW( SCIPcreateExprSum(scip, &expr, 1, &child, &minusone, 0.0, NULL, NULL) );
435 break;
436 }
437
438 case mp::expr::ABS:
439 SCIP_CALL_THROW( SCIPcreateExprAbs(scip, &expr, child, NULL, NULL) );
440 break;
441
442 case mp::expr::POW2:
443 SCIP_CALL_THROW( SCIPcreateExprPow(scip, &expr, child, 2.0, NULL, NULL) );
444 break;
445
446 case mp::expr::SQRT:
447 SCIP_CALL_THROW( SCIPcreateExprPow(scip, &expr, child, 0.5, NULL, NULL) );
448 break;
449
450 case mp::expr::LOG:
451 SCIP_CALL_THROW( SCIPcreateExprLog(scip, &expr, child, NULL, NULL) );
452 break;
453
454 case mp::expr::LOG10: // 1/log(10)*log(child)
455 {
456 SCIP_EXPR* logexpr;
457 SCIP_Real factor = 1.0/log(10.0);
458 SCIP_CALL_THROW( SCIPcreateExprLog(scip, &logexpr, child, NULL, NULL) );
459 SCIP_CALL_THROW( SCIPcreateExprSum(scip, &expr, 1, &logexpr, &factor, 0.0, NULL, NULL) );
460 SCIP_CALL_THROW( SCIPreleaseExpr(scip, &logexpr) );
461 break;
462 }
463
464 case mp::expr::EXP:
465 SCIP_CALL_THROW( SCIPcreateExprExp(scip, &expr, child, NULL, NULL) );
466 break;
467
468 case mp::expr::SIN:
469 SCIP_CALL_THROW( SCIPcreateExprSin(scip, &expr, child, NULL, NULL) );
470 break;
471
472 case mp::expr::COS:
473 SCIP_CALL_THROW( SCIPcreateExprCos(scip, &expr, child, NULL, NULL) );
474 break;
475
476 default:
477 OnUnhandled(mp::expr::str(kind));
478 return NULL;
479 }
480
481 // remember that we have to release this expr
482 exprstorelease.push_back(expr);
483
484 return expr;
485 }
486
487 /// receive notification of a binary expression
488 SCIP_EXPR* OnBinary(
489 mp::expr::Kind kind, ///< expression operand
490 SCIP_EXPR* firstChild, ///< first argument
491 SCIP_EXPR* secondChild ///< second argument
492 )
493 {
494 SCIP_EXPR* expr;
495 SCIP_EXPR* children[2] = { firstChild, secondChild };
496
497 assert(firstChild != NULL);
498 assert(secondChild != NULL);
499
500 switch( kind )
501 {
502 case mp::expr::ADD:
503 SCIP_CALL_THROW( SCIPcreateExprSum(scip, &expr, 2, children, NULL, 0.0, NULL, NULL) );
504 break;
505
506 case mp::expr::SUB:
507 {
508 SCIP_Real coefs[2] = { 1.0, -1.0 };
509 SCIP_CALL_THROW( SCIPcreateExprSum(scip, &expr, 2, children, coefs, 0.0, NULL, NULL) );
510 break;
511 }
512
513 case mp::expr::MUL:
514 SCIP_CALL_THROW( SCIPcreateExprProduct(scip, &expr, 2, children, 1.0, NULL, NULL) );
515 break;
516
517 case mp::expr::DIV:
518 SCIP_CALL_THROW( SCIPcreateExprPow(scip, &children[1], secondChild, -1.0, NULL, NULL) );
519 SCIP_CALL_THROW( SCIPcreateExprProduct(scip, &expr, 2, children, 1.0, NULL, NULL) );
520 SCIP_CALL_THROW( SCIPreleaseExpr(scip, &children[1]) );
521 break;
522
523 case mp::expr::POW_CONST_BASE:
524 case mp::expr::POW_CONST_EXP:
525 case mp::expr::POW:
526 // with some .nl files, we seem to get mp::expr::POW even if base or exponent is constant,
527 // so do not rely on kind but better check expr type
528 if( SCIPisExprValue(scip, secondChild) )
529 {
530 SCIP_CALL_THROW( SCIPcreateExprPow(scip, &expr, firstChild, SCIPgetValueExprValue(secondChild), NULL, NULL) );
531 break;
532 }
533
534 if( SCIPisExprValue(scip, firstChild) && SCIPgetValueExprValue(firstChild) > 0.0 )
535 {
536 // reformulate constant^y as exp(y*log(constant)), if constant > 0.0
537 // if constant < 0, we create an expression and let cons_nonlinear figure out infeasibility somehow
538 SCIP_EXPR* prod;
539
540 SCIP_Real coef = log(SCIPgetValueExprValue(firstChild)); // log(firstChild)
541 SCIP_CALL_THROW( SCIPcreateExprSum(scip, &prod, 1, &secondChild, &coef, 0.0, NULL, NULL) ); // log(firstChild)*secondChild
542 SCIP_CALL_THROW( SCIPcreateExprExp(scip, &expr, prod, NULL, NULL) ); // expr(log(firstChild)*secondChild)
543
544 SCIP_CALL_THROW( SCIPreleaseExpr(scip, &prod) );
545 break;
546 }
547
548 {
549 // reformulate x^y as exp(y*log(x))
550 SCIP_EXPR* prod;
551
552 assert(SCIPisExprValue(scip, secondChild));
553
554 SCIP_CALL_THROW( SCIPcreateExprLog(scip, &children[0], firstChild, NULL, NULL) ); // log(firstChild)
555 SCIP_CALL_THROW( SCIPcreateExprProduct(scip, &prod, 2, children, 1.0, NULL, NULL) ); // log(firstChild)*secondChild
556 SCIP_CALL_THROW( SCIPcreateExprExp(scip, &expr, prod, NULL, NULL) ); // expr(log(firstChild)*secondChild)
557
558 SCIP_CALL_THROW( SCIPreleaseExpr(scip, &prod) );
559 SCIP_CALL_THROW( SCIPreleaseExpr(scip, &children[0]) );
560 break;
561 }
562
563 default:
564 OnUnhandled(mp::expr::str(kind));
565 return NULL;
566 }
567
568 // remember that we have to release this expr
569 exprstorelease.push_back(expr);
570
571 return expr;
572 }
573
574 /// handler to create a list of terms in a sum
575 ///
576 /// NumericArgHandler is copied around, so it keeps only a pointer (with reference counting) to actual data
577 class NumericArgHandler
578 {
579 public:
580 std::shared_ptr<std::vector<SCIP_EXPR*> > v;
581
582 /// constructor
583 NumericArgHandler(
584 int num_args ///< number of terms to expect
585 )
586 : v(new std::vector<SCIP_EXPR*>())
587 {
588 v->reserve(num_args);
589 }
590
591 /// adds term to sum
592 void AddArg(
593 SCIP_EXPR* term ///< term to add
594 )
595 {
596 v->push_back(term);
597 }
598 };
599
600 /// receive notification of the beginning of a summation
601 NumericArgHandler BeginSum(
602 int num_args ///< number of terms to expect
603 )
604 {
605 NumericArgHandler h(num_args);
606 return h;
607 }
608
609 /// receive notification of the end of a summation
610 SCIP_EXPR* EndSum(
611 NumericArgHandler handler ///< handler that handled the sum
612 )
613 {
614 SCIP_EXPR* expr;
615 SCIP_CALL_THROW( SCIPcreateExprSum(scip, &expr, (int)handler.v->size(), handler.v->data(), NULL, 0.0, NULL, NULL) );
616 // remember that we have to release this expr
617 exprstorelease.push_back(expr);
618 return expr;
619 }
620
621 /// receive notification of an objective type and the nonlinear part of an objective expression
622 void OnObj(
623 int objectiveIndex, ///< index of objective
624 mp::obj::Type type, ///< objective sense
625 SCIP_EXPR* nonlinearExpression ///< nonlinear part of objective function
626 )
627 {
628 if( objectiveIndex >= 1 )
629 OnUnhandled("multiple objective functions");
630
631 SCIP_CALL_THROW( SCIPsetObjsense(scip, type == mp::obj::Type::MAX ? SCIP_OBJSENSE_MAXIMIZE : SCIP_OBJSENSE_MINIMIZE) );
632
633 assert(objexpr == NULL);
634
635 if( nonlinearExpression != NULL && SCIPisExprValue(scip, nonlinearExpression) )
636 {
637 // handle objective constant by adding a fixed variable for it
638 SCIP_VAR* objconstvar;
639 SCIP_Real objconst = SCIPgetValueExprValue(nonlinearExpression);
640
641 SCIP_CALL_THROW( SCIPcreateVarBasic(scip, &objconstvar, "objconstant", objconst, objconst, 1.0, SCIP_VARTYPE_CONTINUOUS) );
642 SCIP_CALL_THROW( SCIPaddVar(scip, objconstvar) );
643 SCIP_CALL_THROW( SCIPreleaseVar(scip, &objconstvar) );
644 }
645 else
646 {
647 objexpr = nonlinearExpression;
648 }
649 }
650
651 /// receive notification of an algebraic constraint expression
652 void OnAlgebraicCon(
653 int constraintIndex, ///< index of constraint
654 SCIP_EXPR* expr ///< nonlinear part of constraint
655 )
656 {
657 if( expr != NULL )
658 {
659 SCIP_CALL_THROW( SCIPchgExprNonlinear(scip, probdata->conss[constraintIndex], expr) );
660 }
661 }
662
663 /// handles linear part of a common expression
664 /// sets up a sum expression, if the linear part isn't empty
665 class LinearExprHandler
666 {
667 private:
668 AMPLProblemHandler& amplph;
669 SCIP_EXPR* commonexpr;
670
671 public:
672 /// constructor
673 LinearExprHandler(
674 AMPLProblemHandler& amplph_, ///< problem handler
675 int index, ///< index of common expression
676 int num_linear_terms///< number of terms to expect
677 )
678 : amplph(amplph_),
679 commonexpr(NULL)
680 {
681 if( num_linear_terms > 0 )
682 {
683 SCIP_CALL_THROW( SCIPcreateExprSum(amplph.scip, &commonexpr, 0, NULL, NULL, 0.0, NULL, NULL) );
684 amplph.commonexprs[index] = commonexpr;
685 amplph.exprstorelease.push_back(commonexpr);
686 }
687 }
688
689 /// receives notification of a term in the linear expression
690 void AddTerm(
691 int var_index, ///< AMPL index of variable
692 double coef ///< variable coefficient
693 )
694 {
695 assert(commonexpr != NULL);
696
697 if( coef == 0.0 )
698 return;
699
700 if( var_index < (int)amplph.varexprs.size() )
701 {
702 SCIP_CALL_THROW( SCIPappendExprSumExpr(amplph.scip, commonexpr, amplph.varexprs[var_index], coef) );
703 }
704 else
705 {
706 // the index variable is linear (not sure this can happen here)
707 assert(var_index < amplph.probdata->nvars);
708 SCIP_EXPR* varexpr;
709 SCIP_CALL_THROW( SCIPcreateExprVar(amplph.scip, &varexpr, amplph.probdata->vars[var_index], NULL, NULL) );
710 SCIP_CALL_THROW( SCIPappendExprSumExpr(amplph.scip, commonexpr, varexpr, coef) );
711 SCIP_CALL_THROW( SCIPreleaseExpr(amplph.scip, &varexpr) );
712 }
713 }
714 };
715
716 /// receive notification of the beginning of a common expression (defined variable)
717 LinearExprHandler BeginCommonExpr(
718 int index, ///< index of common expression
719 int num_linear_terms ///< number of terms to expect
720 )
721 {
722 assert(index >= 0);
723 assert(index < (int)commonexprs.size());
724
725 return LinearExprHandler(*this, index, num_linear_terms);
726 }
727
728 /// receive notification of the end of a common expression
729 void EndCommonExpr(
730 int index, ///< index of common expression
731 SCIP_EXPR* expr, ///< nonlinear part of common expression
732 int /* position */ ///< argument that doesn't seem to have any purpose
733 )
734 {
735 if( commonexprs[index] != NULL )
736 {
737 // add expr, if any, to linear part
738 if( expr != NULL )
739 {
740 SCIP_CALL_THROW( SCIPappendExprSumExpr(scip, commonexprs[index], expr, 1.0) );
741 }
742 }
743 else if( expr != NULL )
744 {
745 commonexprs[index] = expr;
746 }
747 }
748
749 /// receive notification of a common expression (defined variable) reference
750 SCIP_EXPR* OnCommonExprRef(
751 int expr_index ///< index of common expression
752 )
753 {
754 assert(expr_index >= 0);
755 assert(expr_index < (int)commonexprs.size());
756 assert(commonexprs[expr_index] != NULL);
757 return commonexprs[expr_index];
758 }
759
760 /// receive notification of variable bounds
761 void OnVarBounds(
762 int variableIndex, ///< AMPL index of variable
763 double variableLB, ///< variable lower bound
764 double variableUB ///< variable upper bound
765 )
766 {
767 assert(variableIndex >= 0);
768 assert(variableIndex < probdata->nvars);
769
770 // as far as I see, ampl::mp gives -inf, +inf for no-bounds, which is always beyond SCIPinfinity()
771 // we ignore bounds outside [-scipinfinity,scipinfinity] here
772 // for binary variables, we also ignore bounds outside [0,1]
773 if( variableLB > (SCIPvarGetType(probdata->vars[variableIndex]) == SCIP_VARTYPE_BINARY ? 0.0 : -SCIPinfinity(scip)) )
774 {
775 SCIP_CALL_THROW( SCIPchgVarLbGlobal(scip, probdata->vars[variableIndex], variableLB) );
776 }
777 if( variableUB < (SCIPvarGetType(probdata->vars[variableIndex]) == SCIP_VARTYPE_BINARY ? 1.0 : SCIPinfinity(scip)) )
778 {
779 SCIP_CALL_THROW( SCIPchgVarUbGlobal(scip, probdata->vars[variableIndex], variableUB) );
780 }
781 }
782
783 /// receive notification of constraint sides
784 void OnConBounds(
785 int index, ///< AMPL index of constraint
786 double lb, ///< constraint left-hand-side
787 double ub ///< constraint right-hand-side
788 )
789 {
790 assert(index >= 0);
791 assert(index < probdata->nconss);
792
793 // nonlinear constraints are first
794 if( index < (int)nlconslin.size() )
795 {
796 if( !SCIPisInfinity(scip, -lb) )
797 {
798 SCIP_CALL_THROW( SCIPchgLhsNonlinear(scip, probdata->conss[index], lb) );
799 }
800 if( !SCIPisInfinity(scip, ub) )
801 {
802 SCIP_CALL_THROW( SCIPchgRhsNonlinear(scip, probdata->conss[index], ub) );
803 }
804 }
805 else
806 {
807 if( !SCIPisInfinity(scip, -lb) )
808 {
809 SCIP_CALL_THROW( SCIPchgLhsLinear(scip, probdata->conss[index], lb) );
810 }
811 if( !SCIPisInfinity(scip, ub) )
812 {
813 SCIP_CALL_THROW( SCIPchgRhsLinear(scip, probdata->conss[index], ub) );
814 }
815 }
816 }
817
818 /// receive notification of the initial value for a variable
819 void OnInitialValue(
820 int var_index, ///< AMPL index of variable
821 double value ///< initial primal value of variable
822 )
823 {
824 if( initsol == NULL )
825 {
826 SCIP_CALL_THROW( SCIPcreateSol(scip, &initsol, NULL) );
827 }
828
829 SCIP_CALL_THROW( SCIPsetSolVal(scip, initsol, probdata->vars[var_index], value) );
830 }
831
832 /// receives notification of the initial value for a dual variable
833 void OnInitialDualValue(
834 int /* con_index */, ///< AMPL index of constraint
835 double /* value */ ///< initial dual value of constraint
836 )
837 {
838 // ignore initial dual value
839 }
840
841 /// receives notification of Jacobian column sizes
842 ColumnSizeHandler OnColumnSizes()
843 {
844 /// use ColumnSizeHandler from upper class, which does nothing
845 return ColumnSizeHandler();
846 }
847
848 /// handling of suffices for variable and constraint flags and SOS constraints
849 ///
850 /// regarding SOS in AMPL, see https://ampl.com/faqs/how-can-i-use-the-solvers-special-ordered-sets-feature/
851 /// we pass the .ref suffix as weight to the SOS constraint handlers
852 /// for a SOS2, the weights determine the order of variables in the set
853 template<typename T> class SuffixHandler
854 {
855 private:
856 AMPLProblemHandler& amplph;
857
858 // type of suffix that is handled, or IGNORE if unsupported suffix
859 enum
860 {
861 IGNORE,
862 CONSINITIAL,
863 CONSSEPARATE,
864 CONSENFORCE,
865 CONSCHECK,
866 CONSPROPAGATE,
867 CONSDYNAMIC,
868 CONSREMOVABLE,
869 VARINITIAL,
870 VARREMOVABLE,
871 VARSOSNO,
872 VARREF,
873 } suffix;
874
875 public:
876 /// constructor
877 SuffixHandler(
878 AMPLProblemHandler& amplph_, ///< problem handler
879 fmt::StringRef name, ///< name of suffix
880 mp::suf::Kind kind ///< whether suffix applies to var, cons, etc
881 )
882 : amplph(amplph_),
883 suffix(IGNORE)
884 {
885 switch( kind )
886 {
887 case mp::suf::Kind::CON:
888 if( strncmp(name.data(), "initial", name.size()) == 0 )
889 {
890 suffix = CONSINITIAL;
891 }
892 else if( strncmp(name.data(), "separate", name.size()) == 0 )
893 {
894 suffix = CONSSEPARATE;
895 }
896 else if( strncmp(name.data(), "enforce", name.size()) == 0 )
897 {
898 suffix = CONSENFORCE;
899 }
900 else if( strncmp(name.data(), "check", name.size()) == 0 )
901 {
902 suffix = CONSCHECK;
903 }
904 else if( strncmp(name.data(), "propagate", name.size()) == 0 )
905 {
906 suffix = CONSPROPAGATE;
907 }
908 else if( strncmp(name.data(), "dynamic", name.size()) == 0 )
909 {
910 suffix = CONSDYNAMIC;
911 }
912 else if( strncmp(name.data(), "removable", name.size()) == 0 )
913 {
914 suffix = CONSREMOVABLE;
915 }
916 else
917 {
918 SCIPverbMessage(amplph.scip, SCIP_VERBLEVEL_HIGH, NULL, "Unknown constraint suffix <%.*s>. Ignoring.\n", (int)name.size(), name.data());
919 }
920 break;
921
922 case mp::suf::Kind::VAR:
923 {
924 if( strncmp(name.data(), "initial", name.size()) == 0 )
925 {
926 suffix = VARINITIAL;
927 }
928 else if( strncmp(name.data(), "removable", name.size()) == 0 )
929 {
930 suffix = VARREMOVABLE;
931 }
932 else if( strncmp(name.data(), "sosno", name.size()) == 0 )
933 {
934 // SOS membership
935 suffix = VARSOSNO;
936 }
937 else if( strncmp(name.data(), "ref", name.size()) == 0 )
938 {
939 // SOS weights
940 suffix = VARREF;
941 amplph.sosweights.resize(amplph.probdata->nvars, 0);
942 }
943 else
944 {
945 SCIPverbMessage(amplph.scip, SCIP_VERBLEVEL_HIGH, NULL, "Unknown variable suffix <%.*s>. Ignoring.\n", (int)name.size(), name.data());
946 }
947 break;
948
949 case mp::suf::Kind::OBJ:
950 SCIPverbMessage(amplph.scip, SCIP_VERBLEVEL_HIGH, NULL, "Unknown objective suffix <%.*s>. Ignoring.\n", (int)name.size(), name.data());
951 break;
952
953 case mp::suf::Kind::PROBLEM:
954 SCIPverbMessage(amplph.scip, SCIP_VERBLEVEL_HIGH, NULL, "Unknown problem suffix <%.*s>. Ignoring.\n", (int)name.size(), name.data());
955 break;
956 }
957 }
958 }
959
960 void SetValue(
961 int index, ///< index of variable, constraint, etc
962 T value ///< value of suffix
963 )
964 {
965 assert(index >= 0);
966 switch( suffix )
967 {
968 case IGNORE :
969 return;
970
971 case CONSINITIAL:
972 SCIP_CALL_THROW( SCIPsetConsInitial(amplph.scip, amplph.probdata->conss[index], value == 1) );
973 break;
974
975 case CONSSEPARATE:
976 SCIP_CALL_THROW( SCIPsetConsSeparated(amplph.scip, amplph.probdata->conss[index], value == 1) );
977 break;
978
979 case CONSENFORCE:
980 SCIP_CALL_THROW( SCIPsetConsEnforced(amplph.scip, amplph.probdata->conss[index], value == 1) );
981 break;
982
983 case CONSCHECK:
984 SCIP_CALL_THROW( SCIPsetConsChecked(amplph.scip, amplph.probdata->conss[index], value == 1) );
985 break;
986
987 case CONSPROPAGATE:
988 SCIP_CALL_THROW( SCIPsetConsPropagated(amplph.scip, amplph.probdata->conss[index], value == 1) );
989 break;
990
991 case CONSDYNAMIC:
992 SCIP_CALL_THROW( SCIPsetConsDynamic(amplph.scip, amplph.probdata->conss[index], value == 1) );
993 break;
994
995 case CONSREMOVABLE:
996 SCIP_CALL_THROW( SCIPsetConsRemovable(amplph.scip, amplph.probdata->conss[index], value == 1) );
997 break;
998
999 case VARINITIAL:
1000 assert(index < amplph.probdata->nvars);
1001 SCIP_CALL_THROW( SCIPvarSetInitial(amplph.probdata->vars[index], value == 1) );
1002 break;
1003
1004 case VARREMOVABLE:
1005 assert(index < amplph.probdata->nvars);
1006 SCIP_CALL_THROW( SCIPvarSetRemovable(amplph.probdata->vars[index], value == 1) );
1007 break;
1008
1009 case VARSOSNO:
1010 // remember that variable index belongs to SOS identified by value
1011 amplph.sosvars[(int)value].push_back(index);
1012 break;
1013
1014 case VARREF:
1015 // remember that variable index has weight value
1016 amplph.sosweights[index] = (int)value;
1017 break;
1018 }
1019 }
1020 };
1021
1022 typedef SuffixHandler<int> IntSuffixHandler;
1023 /// receive notification of an integer suffix
1024 IntSuffixHandler OnIntSuffix(
1025 fmt::StringRef name, ///< suffix name, not null-terminated
1026 mp::suf::Kind kind, ///< suffix kind
1027 int /*num_values*/ ///< number of values to expect
1028 )
1029 {
1030 return IntSuffixHandler(*this, name, kind);
1031 }
1032
1033 typedef SuffixHandler<SCIP_Real> DblSuffixHandler;
1034 /// receive notification of a double suffix
1035 DblSuffixHandler OnDblSuffix(
1036 fmt::StringRef name, ///< suffix name, not null-terminated
1037 mp::suf::Kind kind, ///< suffix kind
1038 int /*num_values*/ ///< number of values to expect
1039 )
1040 {
1041 return DblSuffixHandler(*this, name, kind);
1042 }
1043
1044 /// handles receiving the linear part of an objective or constraint
1045 ///
1046 /// for objective, set the objective-coefficient of the variable
1047 /// for linear constraints, add to the constraint
1048 /// for nonlinear constraints, add to nlconslin vector; adding to constraint later
1049 class LinearPartHandler
1050 {
1051 private:
1052 AMPLProblemHandler& amplph;
1053 int constraintIndex;
1054
1055 public:
1056 // constructor for constraint
1057 explicit LinearPartHandler(
1058 AMPLProblemHandler& amplph_, ///< problem handler
1059 int constraintIndex_///< constraint index
1060 )
1061 : amplph(amplph_),
1062 constraintIndex(constraintIndex_)
1063 {
1064 assert(constraintIndex_ >= 0);
1065 assert(constraintIndex_ < amplph.probdata->nconss);
1066 }
1067
1068 // constructor for linear objective
1069 explicit LinearPartHandler(
1070 AMPLProblemHandler& amplph_ ///< problem handler
1071 )
1072 : amplph(amplph_),
1073 constraintIndex(-1)
1074 { }
1075
1076 void AddTerm(
1077 int variableIndex, ///< AMPL index of variable
1078 double coefficient ///< coefficient of variable
1079 )
1080 {
1081 assert(variableIndex >= 0);
1082 assert(variableIndex < amplph.probdata->nvars);
1083
1084 if( coefficient == 0.0 )
1085 return;
1086
1087 if( constraintIndex < 0 )
1088 {
1089 SCIP_CALL_THROW( SCIPchgVarObj(amplph.scip, amplph.probdata->vars[variableIndex], coefficient) );
1090 }
1091 else if( constraintIndex < (int)amplph.nlconslin.size() )
1092 {
1093 amplph.nlconslin[constraintIndex].push_back(std::pair<SCIP_Real, SCIP_VAR*>(coefficient, amplph.probdata->vars[variableIndex]));
1094 }
1095 else
1096 {
1097 SCIP_CONS* lincons = amplph.probdata->conss[constraintIndex];
1098 SCIP_CALL_THROW( SCIPaddCoefLinear(amplph.scip, lincons, amplph.probdata->vars[variableIndex], coefficient) );
1099 }
1100 }
1101 };
1102
1103 typedef LinearPartHandler LinearObjHandler;
1104
1105 /// receive notification of the linear part of an objective
1106 LinearPartHandler OnLinearObjExpr(
1107 int objectiveIndex, ///< index of objective
1108 int /* numLinearTerms *////< number of terms to expect
1109 )
1110 {
1111 if( objectiveIndex >= 1 )
1112 OnUnhandled("multiple objective functions");
1113
1114 return LinearObjHandler(*this);
1115 }
1116
1117 typedef LinearPartHandler LinearConHandler;
1118
1119 /// receive notification of the linear part of a constraint
1120 LinearConHandler OnLinearConExpr(
1121 int constraintIndex, ///< index of constraint
1122 int /* numLinearTerms *////< number of terms to expect
1123 )
1124 {
1125 return LinearConHandler(*this, constraintIndex);
1126 }
1127
1128 /// receive notification of the end of the input
1129 ///
1130 /// - setup all nonlinear constraints and add them to SCIP
1131 /// - add linear constraints to SCIP (should be after nonlinear ones to respect order in .nl file)
1132 /// - add initial solution, if initial values were given
1133 void EndInput()
1134 {
1135 // turn nonlinear objective into constraint
1136 // min f(x) -> min z s.t. f(x) - z <= 0
1137 // max f(x) -> max z s.t. 0 <= f(x) - z
1138 if( objexpr != NULL )
1139 {
1140 SCIP_CONS* objcons;
1141 SCIP_VAR* objvar;
1142
1143 SCIP_CALL_THROW( SCIPcreateVarBasic(scip, &objvar, "nlobjvar", -SCIPinfinity(scip), SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) );
1144 SCIP_CALL_THROW( SCIPaddVar(scip, objvar) );
1145
1146 SCIP_CALL_THROW( SCIPcreateConsBasicNonlinear(scip, &objcons, "objcons", objexpr,
1147 SCIPgetObjsense(scip) == SCIP_OBJSENSE_MINIMIZE ? -SCIPinfinity(scip) : 0.0,
1148 SCIPgetObjsense(scip) == SCIP_OBJSENSE_MAXIMIZE ? SCIPinfinity(scip) : 0.0) );
1149 SCIP_CALL_THROW( SCIPaddLinearVarNonlinear(scip, objcons, objvar, -1.0) );
1150 SCIP_CALL_THROW( SCIPaddCons(scip, objcons) );
1151
1152 SCIP_CALL_THROW( SCIPreleaseCons(scip, &objcons) );
1153 SCIP_CALL_THROW( SCIPreleaseVar(scip, &objvar) );
1154 }
1155
1156 // add linear terms to expressions of nonlinear constraints (should be ok to do this one-by-one for now)
1157 for( size_t i = 0; i < nlconslin.size(); ++i )
1158 {
1159 for( size_t j = 0; j < nlconslin[i].size(); ++j )
1160 {
1161 SCIP_CALL_THROW( SCIPaddLinearVarNonlinear(scip, probdata->conss[i], nlconslin[i][j].second, nlconslin[i][j].first) );
1162 }
1163 }
1164
1165 // add constraints
1166 for( int i = 0; i < probdata->nconss; ++i )
1167 {
1168 SCIP_CALL_THROW( SCIPaddCons(scip, probdata->conss[i]) );
1169 }
1170
1171 // add SOS constraints
1172 std::vector<SCIP_VAR*> setvars; // variables in one SOS
1173 std::vector<SCIP_Real> setweights; // weights for one SOS
1174 if( !sosvars.empty() )
1175 {
1176 setvars.resize(probdata->nvars);
1177 probdata->islp = false;
1178 }
1179 if( !sosweights.empty() )
1180 setweights.resize(probdata->nvars);
1181 for( std::map<int, std::vector<int> >::iterator sosit(sosvars.begin()); sosit != sosvars.end(); ++sosit )
1182 {
1183 assert(sosit->first != 0);
1184 assert(!sosit->second.empty());
1185
1186 // a negative SOS identifier means SOS2
1187 bool issos2 = sosit->first < 0;
1188
1189 if( issos2 && sosweights.empty() )
1190 {
1191 // if no .ref suffix was given for a SOS2 constraint, then we consider this as an error
1192 // since the weights determine the order
1193 // for a SOS1, the weights only specify branching preference, so can treat them as optional
1194 OnUnhandled("SOS2 requires variable .ref suffix");
1195 }
1196
1197 for( size_t i = 0; i < sosit->second.size(); ++i )
1198 {
1199 int varidx = sosit->second[i];
1200 setvars[i] = probdata->vars[varidx];
1201
1202 if( issos2 && sosweights[varidx] == 0 )
1203 // 0 is the default if no ref was given for a variable; we don't allow this for SOS2
1204 OnUnhandled("Missing .ref value for SOS2 variable");
1205 if( !sosweights.empty() )
1206 setweights[i] = (SCIP_Real)sosweights[varidx];
1207 }
1208
1209 SCIP_CONS* cons;
1210 char name[20];
1211 if( !issos2 )
1212 {
|
(1) Event secure_coding: |
[VERY RISKY]. Using "sprintf(char * restrict, char const * restrict, ...)" can cause a buffer overflow when done incorrectly. Because sprintf() assumes an arbitrarily long string, callers must be careful not to overflow the actual space of the destination. Use snprintf() instead, or correct precision specifiers. |
| Also see events: |
[secure_coding] |
1213 (void) sprintf(name, "sos1_%d", sosit->first);
1214 SCIP_CALL_THROW( SCIPcreateConsBasicSOS1(scip, &cons, name, sosit->second.size(), setvars.data(), setweights.empty() ? NULL : setweights.data()) );
1215 }
1216 else
1217 {
|
(2) Event secure_coding: |
[VERY RISKY]. Using "sprintf(char * restrict, char const * restrict, ...)" can cause a buffer overflow when done incorrectly. Because sprintf() assumes an arbitrarily long string, callers must be careful not to overflow the actual space of the destination. Use snprintf() instead, or correct precision specifiers. |
| Also see events: |
[secure_coding] |
1218 (void) sprintf(name, "sos2_%d", -sosit->first);
1219 SCIP_CALL_THROW( SCIPcreateConsBasicSOS2(scip, &cons, name, sosit->second.size(), setvars.data(), setweights.data()) );
1220 }
1221 SCIP_CALL_THROW( SCIPaddCons(scip, cons) );
1222 SCIP_CALL_THROW( SCIPreleaseCons(scip, &cons) );
1223 }
1224
1225 // add initial solution
1226 if( initsol != NULL )
1227 {
1228 SCIP_Bool stored;
1229 SCIP_CALL_THROW( SCIPaddSolFree(scip, &initsol, &stored) );
1230 }
1231
1232 // release expressions
1233 SCIP_CALL_THROW( cleanup() );
1234 }
1235
1236 /// releases expressions and linear constraints from data
1237 ///
1238 /// should be called if there was an error while reading the .nl file
1239 /// this is not in the destructor, because we want to return SCIP_RETCODE
1240 SCIP_RETCODE cleanup()
1241 {
1242 // release initial sol (in case EndInput() wasn't called)
1243 if( initsol != NULL )
1244 {
1245 SCIP_CALL( SCIPfreeSol(scip, &initsol) );
1246 }
1247
1248 // release created expressions (they should all be used in other expressions or constraints now)
1249 while( !exprstorelease.empty() )
1250 {
1251 SCIP_CALL( SCIPreleaseExpr(scip, &exprstorelease.back()) );
1252 exprstorelease.pop_back();
1253 }
1254
1255 // release variable expressions (they should all be used in other expressions or constraints now)
1256 while( !varexprs.empty() )
1257 {
1258 SCIP_CALL( SCIPreleaseExpr(scip, &varexprs.back()) );
1259 varexprs.pop_back();
1260 }
1261
1262 return SCIP_OKAY;
1263 }
1264 };
1265
1266
1267 /*
1268 * Callback methods of probdata
1269 */
1270
1271 /** frees user data of original problem (called when the original problem is freed) */
1272 static
1273 SCIP_DECL_PROBDELORIG(probdataDelOrigNl)
1274 {
1275 int i;
1276
1277 assert((*probdata)->vars != NULL || (*probdata)->nvars == 0);
1278 assert((*probdata)->conss != NULL || (*probdata)->conss == 0);
1279
1280 for( i = 0; i < (*probdata)->nconss; ++i )
1281 {
1282 SCIP_CALL( SCIPreleaseCons(scip, &(*probdata)->conss[i]) );
1283 }
1284 SCIPfreeBlockMemoryArrayNull(scip, &(*probdata)->conss, (*probdata)->nconss);
1285
1286 for( i = 0; i < (*probdata)->nvars; ++i )
1287 {
1288 SCIP_CALL( SCIPreleaseVar(scip, &(*probdata)->vars[i]) );
1289 }
1290 SCIPfreeBlockMemoryArrayNull(scip, &(*probdata)->vars, (*probdata)->nvars);
1291
1292 SCIPfreeBlockMemoryArrayNull(scip, &(*probdata)->filenamestub, (*probdata)->filenamestublen+5);
1293
1294 SCIPfreeMemory(scip, probdata);
1295
1296 return SCIP_OKAY;
1297 }
1298
1299 /*
1300 * Callback methods of reader
1301 */
1302
1303 /** copy method for reader plugins (called when SCIP copies plugins) */
1304 static
1305 SCIP_DECL_READERCOPY(readerCopyNl)
1306 { /*lint --e{715}*/
1307 assert(scip != NULL);
1308
1309 SCIP_CALL( SCIPincludeReaderNl(scip) );
1310
1311 return SCIP_OKAY;
1312 }
1313
1314 /** problem reading method of reader */
1315 static
1316 SCIP_DECL_READERREAD(readerReadNl)
1317 { /*lint --e{715}*/
1318 assert(scip != NULL);
1319 assert(reader != NULL);
1320 assert(filename != NULL);
1321 assert(result != NULL);
1322
1323 try
1324 {
1325 // try to read the .nl file and setup SCIP problem
1326 AMPLProblemHandler handler(scip, filename);
1327 try
1328 {
1329 mp::ReadNLFile(filename, handler);
1330 }
1331 catch( const mp::UnsupportedError& e )
1332 {
1333 SCIPerrorMessage("unsupported construct in AMPL .nl file %s: %s\n", filename, e.what());
1334
1335 SCIP_CALL( handler.cleanup() );
1336
1337 return SCIP_READERROR;
1338 }
1339 catch( const mp::Error& e )
1340 {
1341 // some other error from ampl/mp, maybe invalid .nl file
1342 SCIPerrorMessage("%s\n", e.what());
1343
1344 SCIP_CALL( handler.cleanup() );
1345
1346 return SCIP_READERROR;
1347 }
1348 catch( const fmt::SystemError& e )
1349 {
1350 // probably a file open error, probably because file not found
1351 SCIPerrorMessage("%s\n", e.what());
1352
1353 SCIP_CALL( handler.cleanup() );
1354
1355 return SCIP_NOFILE;
1356 }
1357 catch( const std::bad_alloc& e )
1358 {
1359 SCIPerrorMessage("Out of memory: %s\n", e.what());
1360
1361 SCIP_CALL( handler.cleanup() );
1362
1363 return SCIP_NOMEMORY;
1364 }
1365 }
1366 catch( const std::exception& e )
1367 {
1368 SCIPerrorMessage("%s\n", e.what());
1369 return SCIP_ERROR;
1370 }
1371
1372 *result = SCIP_SUCCESS;
1373
1374 return SCIP_OKAY;
1375 }
1376
1377 /*
1378 * reader specific interface methods
1379 */
1380
1381 /** includes the AMPL .nl file reader in SCIP */
1382 SCIP_RETCODE SCIPincludeReaderNl(
1383 SCIP* scip /**< SCIP data structure */
1384 )
1385 {
1386 SCIP_READER* reader = NULL;
1387
1388 /* include reader */
1389 SCIP_CALL( SCIPincludeReaderBasic(scip, &reader, READER_NAME, READER_DESC, READER_EXTENSION, NULL) );
1390 assert(reader != NULL);
1391
1392 /* set non fundamental callbacks via setter functions */
1393 SCIP_CALL( SCIPsetReaderCopy(scip, reader, readerCopyNl) );
1394 SCIP_CALL( SCIPsetReaderRead(scip, reader, readerReadNl) );
1395
1396 SCIP_CALL( SCIPincludeExternalCodeInformation(scip, "AMPL/MP 4e2d45c4", "AMPL .nl file reader library (github.com/ampl/mp)") );
1397
1398 return SCIP_OKAY;
1399 }
1400
1401 /** writes AMPL solution file
1402 *
1403 * problem must have been read with .nl reader
1404 */
1405 SCIP_RETCODE SCIPwriteSolutionNl(
1406 SCIP* scip /**< SCIP data structure */
1407 )
1408 {
1409 SCIP_PROBDATA* probdata;
1410
1411 assert(scip != NULL);
1412
1413 probdata = SCIPgetProbData(scip);
1414 if( probdata == NULL )
1415 {
1416 SCIPerrorMessage("No AMPL nl file read. Cannot write AMPL solution.\n");
1417 return SCIP_ERROR;
1418 }
1419
1420 probdata->filenamestub[probdata->filenamestublen] = '.';
1421 probdata->filenamestub[probdata->filenamestublen+1] = 's';
1422 probdata->filenamestub[probdata->filenamestublen+2] = 'o';
1423 probdata->filenamestub[probdata->filenamestublen+3] = 'l';
1424 probdata->filenamestub[probdata->filenamestublen+4] = '\0';
1425
1426 FILE* solfile = fopen(probdata->filenamestub, "w");
1427 if( solfile == NULL )
1428 {
1429 SCIPerrorMessage("could not open file <%s> for writing\n", probdata->filenamestub);
1430 probdata->filenamestub[probdata->filenamestublen] = '\0';
1431
1432 return SCIP_WRITEERROR;
1433 }
1434 probdata->filenamestub[probdata->filenamestublen] = '\0';
1435
1436 // see ampl/mp:sol.h:WriteSolFile() (seems buggy, https://github.com/ampl/mp/issues/135) and asl/writesol.c for solution file format
1437 SCIP_CALL( SCIPprintStatus(scip, solfile) );
1438 SCIPinfoMessage(scip, solfile, "\n\n");
1439
1440 SCIPinfoMessage(scip, solfile, "Options\n%d\n", probdata->namplopts);
1441 for( int i = 0; i < probdata->namplopts; ++i )
1442 SCIPinfoMessage(scip, solfile, "%d\n", probdata->amplopts[i]);
1443
1444 bool haveprimal = SCIPgetBestSol(scip) != NULL;
1445 bool havedual = probdata->islp && SCIPgetStage(scip) == SCIP_STAGE_SOLVED && !SCIPhasPerformedPresolve(scip);
1446
1447 SCIPinfoMessage(scip, solfile, "%d\n%d\n", probdata->nconss, havedual ? probdata->nconss : 0);
1448 SCIPinfoMessage(scip, solfile, "%d\n%d\n", probdata->nvars, haveprimal ? probdata->nvars : 0);
1449
1450 SCIPdebug( SCIPprintSol(scip, SCIPgetBestSol(scip), NULL, TRUE); )
1451
1452 if( havedual )
1453 for( int c = 0; c < probdata->nconss; ++c )
1454 {
1455 SCIP_CONS* transcons;
1456 SCIP_Real dualval;
1457
1458 /* dual solution is created by LP solver and therefore only available for linear constraints */
1459 SCIP_CALL( SCIPgetTransformedCons(scip, probdata->conss[c], &transcons) );
1460 assert(transcons == NULL || strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(transcons)), "linear") == 0);
1461
1462 if( transcons == NULL )
1463 dualval = 0.0;
1464 else if( SCIPgetObjsense(scip) == SCIP_OBJSENSE_MINIMIZE )
1465 dualval = SCIPgetDualsolLinear(scip, transcons);
1466 else
1467 dualval = -SCIPgetDualsolLinear(scip, transcons);
1468 assert(dualval != SCIP_INVALID);
1469
1470 SCIPinfoMessage(scip, solfile, "%.17g\n", dualval);
1471 }
1472
1473 if( haveprimal )
1474 for( int i = 0; i < probdata->nvars; ++i )
1475 SCIPinfoMessage(scip, solfile, "%.17g\n", SCIPgetSolVal(scip, SCIPgetBestSol(scip), probdata->vars[i]));
1476
1477 /* AMPL solve status codes are at http://www.ampl.com/NEW/statuses.html
1478 * number string interpretation
1479 * 0 - 99 solved optimal solution found
1480 * 100 - 199 solved? optimal solution indicated, but error likely
1481 * 200 - 299 infeasible constraints cannot be satisfied
1482 * 300 - 399 unbounded objective can be improved without limit
1483 * 400 - 499 limit stopped by a limit that you set (such as on iterations)
1484 * 500 - 599 failure stopped by an error condition in the solver routines
1485 */
1486 int solve_result_num;
1487 switch( SCIPgetStatus(scip) )
1488 {
1489 case SCIP_STATUS_UNKNOWN:
1490 solve_result_num = 500;
1491 break;
1492 case SCIP_STATUS_USERINTERRUPT:
1493 solve_result_num = 450;
1494 break;
1495 case SCIP_STATUS_NODELIMIT:
1496 solve_result_num = 400;
1497 break;
1498 case SCIP_STATUS_TOTALNODELIMIT:
1499 solve_result_num = 401;
1500 break;
1501 case SCIP_STATUS_STALLNODELIMIT:
1502 solve_result_num = 402;
1503 break;
1504 case SCIP_STATUS_TIMELIMIT:
1505 solve_result_num = 403;
1506 break;
1507 case SCIP_STATUS_MEMLIMIT:
1508 solve_result_num = 404;
1509 break;
1510 case SCIP_STATUS_GAPLIMIT:
1511 solve_result_num = 405;
1512 break;
1513 case SCIP_STATUS_SOLLIMIT:
1514 solve_result_num = 406;
1515 break;
1516 case SCIP_STATUS_BESTSOLLIMIT:
1517 solve_result_num = 407;
1518 break;
1519 case SCIP_STATUS_OPTIMAL:
1520 solve_result_num = 0;
1521 break;
1522 case SCIP_STATUS_INFEASIBLE:
1523 solve_result_num = 200;
1524 break;
1525 case SCIP_STATUS_UNBOUNDED:
1526 solve_result_num = 300;
1527 break;
1528 case SCIP_STATUS_INFORUNBD:
1529 solve_result_num = 299;
1530 break;
1531 default:
1532 /* solve_result_num = 500; */
1533 SCIPerrorMessage("invalid status code <%d>\n", SCIPgetStatus(scip));
1534 return SCIP_INVALIDDATA;
1535 }
1536 SCIPinfoMessage(scip, solfile, "objno 0 %d\n", solve_result_num);
1537
1538 if( fclose(solfile) != 0 )
1539 {
1540 SCIPerrorMessage("could not close solution file after writing\n");
1541 return SCIP_WRITEERROR;
1542 }
1543
1544 return SCIP_OKAY;
1545 }
1546