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