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