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
|
(1) Event missing_assign: |
Class "AMPLProblemHandler" owns resources that are freed in its destructor but has no user-written assignment operator. |
|
(2) Event free_resource: |
The destructor frees member "colfile". [details] |
|
(3) Event free_resource: |
The destructor frees member "rowfile". [details] |
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
|
(1) Event freed_arg: |
"operator delete" frees parameter "this->colfile". |
257 delete colfile;
|
(1) Event freed_arg: |
"operator delete" frees parameter "this->rowfile". |
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 {
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 {
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