diff --git a/CHANGELOG b/CHANGELOG index 96810bdd63..959d23454d 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -30,10 +30,13 @@ Build system ### Makefile +- the check/solchecker build now finds a GMP installation by Homebrew on macOS/arm64 + Miscellaneous ------------- -- The output precision for writing CIP/POB files has been increased at several places for nonlinear constraints. +- the output precision for writing CIP/OPB files has been increased at several places for nonlinear constraints +- the solchecker tool now also supports SCIP solution format Known bugs ---------- diff --git a/check/solchecker/Makefile b/check/solchecker/Makefile index 3bc26bfbe4..8c38b3860b 100644 --- a/check/solchecker/Makefile +++ b/check/solchecker/Makefile @@ -14,8 +14,8 @@ LINKCXX_l = -l LINKCXX_o = -o # the trailing space is important FLAGS = OFLAGS = -CXXFLAGS = -I/usr/local/include -LDFLAGS = -L/usr/local/lib +CXXFLAGS = -I/usr/local/include -I/opt/homebrew/include +LDFLAGS = -L/usr/local/lib -L/opt/homebrew/lib DOXY = doxygen diff --git a/check/solchecker/README b/check/solchecker/README index 5a1539c8f6..6d04cb07a8 100644 --- a/check/solchecker/README +++ b/check/solchecker/README @@ -3,7 +3,7 @@ Solution Checker This code implements a solution checker for MIP models. It reads a MIP model from a (gzipped) MPS file and a -solution from a SOL file and checks the feasibility of this +solution from a non-gzipped SOL file and checks the feasibility of this solution. All input is internally stored as arbitrary precision numbers (thanks to GNU GMP) and all calculations are done in arbitrary precision. The code checks the integrality, the diff --git a/check/solchecker/src/gmputils.cpp b/check/solchecker/src/gmputils.cpp index 6f4e0dbe5d..0780f785e2 100644 --- a/check/solchecker/src/gmputils.cpp +++ b/check/solchecker/src/gmputils.cpp @@ -1,6 +1,6 @@ /** * @file gmputils.cpp - * @brief Basic classes to for rational arithmetic + * @brief Basic classes for rational arithmetic * * @author Domenico Salvagnin * @author Thorsten Koch @@ -11,6 +11,7 @@ #include #include + char Rational::buffer[] = {'\0'}; Rational::Rational(int num, int den) @@ -40,7 +41,7 @@ Rational::~Rational() Rational& Rational::operator=(const Rational& rhs) { - if (this != &rhs) + if( this != &rhs ) mpq_set(number, rhs.number); return *this; } @@ -116,7 +117,7 @@ void Rational::integralityViolation(Rational& violation) const } // otherwise, we must check w.r.t. the given tolerance // first calculate the fractional part - violation = (*this); + violation = *this; violation.abs(); mpz_t r; mpz_init(r); @@ -136,7 +137,8 @@ void Rational::toZero() bool Rational::isInteger(const Rational& tolerance) const { // if denominator is 1, then it is an integer for sure - if (mpz_cmp_ui(mpq_denref(number), 1) == 0) return true; + if( mpz_cmp_ui(mpq_denref(number), 1) == 0 ) + return true; // otherwise, we must check w.r.t. the given tolerance // first calculate the fractional part Rational viol(*this); @@ -170,39 +172,38 @@ bool Rational::isZero() const void Rational::fromString(const char* num) { char* tmp = &buffer[0]; - int k = 0; - int exponent = 0; - int fraction = 0; + int exponent = 0; + int fraction = 0; + int k = 0; assert(num != NULL); - assert(strlen(num) < 32); - // Skip initial whitespace - while(isspace(*num)) + // skip initial whitespaces + while( isspace(*num) ) num++; - // Skip initial +/- - if (*num == '+') + // skip initial sign + if( *num == '+' ) num++; - else if (*num == '-') + else if( *num == '-' ) tmp[k++] = *num++; - for(int i = 0; num[i] != '\0'; i++) + for( int i = 0; num[i] != '\0'; ++i ) { - if (isdigit(num[i])) + if( isdigit(num[i]) ) { tmp[k++] = num[i]; exponent -= fraction; } - else if (num[i] == '.') + else if( num[i] == '.' ) fraction = 1; - else if (tolower(num[i]) == 'e') + else if( tolower(num[i]) == 'e' ) { exponent += atoi(&num[i + 1]); break; } } - while(exponent > 0) + while( exponent > 0 ) { tmp[k++] = '0'; exponent--; @@ -210,7 +211,7 @@ void Rational::fromString(const char* num) tmp[k++] = '/'; tmp[k++] = '1'; - while(exponent < 0) + while( exponent < 0 ) { tmp[k++] = '0'; exponent++; diff --git a/check/solchecker/src/gmputils.h b/check/solchecker/src/gmputils.h index 9d22c697fb..764bd155bd 100644 --- a/check/solchecker/src/gmputils.h +++ b/check/solchecker/src/gmputils.h @@ -1,6 +1,6 @@ /** * @file gmputils.h - * @brief Basic classes to for rational arithmetic + * @brief Basic classes for rational arithmetic * * @author Domenico Salvagnin */ @@ -11,6 +11,7 @@ #include #include + /** * @brief Simple wrapper class around GMP mpq_t. * Not a full-blown wrapper such as the ones in GMP++, but should be enough for the job. diff --git a/check/solchecker/src/model.cpp b/check/solchecker/src/model.cpp index dd1409ef36..d3e3546f9c 100644 --- a/check/solchecker/src/model.cpp +++ b/check/solchecker/src/model.cpp @@ -15,6 +15,7 @@ #define SOL_MAX_LINELEN 1024 #define BLANK ' ' + Var::Var(const char* _name, VarType _type, const Rational& _lb, const Rational& _ub, const Rational& _obj) : name(_name), type(_type), lb(_lb), ub(_ub), objCoef(_obj) {} @@ -41,7 +42,7 @@ bool Var::checkBounds(const Rational& boundTolerance) const relaxedLb -= lbtol; Rational relaxedUb(ub); relaxedUb += ubtol; - if( (value < relaxedLb) || (value > relaxedUb) ) + if( value < relaxedLb || value > relaxedUb ) { std::cerr << "Failed bound check: "; print(std::cerr); @@ -53,7 +54,7 @@ bool Var::checkBounds(const Rational& boundTolerance) const bool Var::checkIntegrality(const Rational& intTolerance) const { - if( (type != CONTINUOUS) && !value.isInteger(intTolerance) ) + if( type != CONTINUOUS && !value.isInteger(intTolerance) ) { std::cerr << "Failed integrality check: "; print(std::cerr); @@ -68,16 +69,19 @@ void Var::boundsViolation(Rational& boundViol) const Rational lbViol; Rational ubViol; sub(lbViol, lb, value); - if( lbViol.isNegative() ) lbViol.toZero(); + if( lbViol.isNegative() ) + lbViol.toZero(); sub(ubViol, value, ub); - if( ubViol.isNegative() ) ubViol.toZero(); + if( ubViol.isNegative() ) + ubViol.toZero(); max(boundViol, lbViol, ubViol); } void Var::integralityViolation(Rational& intViol) const { intViol.toZero(); - if( type == CONTINUOUS ) return; + if( type == CONTINUOUS ) + return; value.integralityViolation(intViol); } @@ -125,8 +129,10 @@ bool LinearConstraint::check(const Rational& tolerance) const { Rational prod; mult(prod, coefs[i], vars[i]->value); - if (prod.isPositive()) posact += prod; - else negact += prod; + if( prod.isPositive() ) + posact += prod; + else + negact += prod; } Rational activity; activity += posact; @@ -154,7 +160,7 @@ bool LinearConstraint::check(const Rational& tolerance) const Rational relaxedRhs(rhs); relaxedRhs += rhstol; // check lhs and rhs - if (activity < relaxedLhs || activity > relaxedRhs) + if( activity < relaxedLhs || activity > relaxedRhs ) { std::cerr << std::setprecision(16) << "Failed check for cons " << name << ": " << activity.toDouble() << " not in [" << relaxedLhs.toDouble() << "," << relaxedRhs.toDouble() << "] -- Exact wrt linear tol: " @@ -175,9 +181,11 @@ void LinearConstraint::violation(Rational& viol) const Rational lhsViol; Rational rhsViol; sub(lhsViol, lhs, activity); - if( lhsViol.isNegative() ) lhsViol.toZero(); + if( lhsViol.isNegative() ) + lhsViol.toZero(); sub(rhsViol, activity, rhs); - if( rhsViol.isNegative() ) rhsViol.toZero(); + if( rhsViol.isNegative() ) + rhsViol.toZero(); max(viol, lhsViol, rhsViol); } @@ -185,9 +193,7 @@ void LinearConstraint::print(std::ostream& out) const { out << name << " " << type << ": " << lhs.toString() << " <= "; for( unsigned int i = 0; i < vars.size(); ++i ) - { out << coefs[i].toString() << " " << vars[i]->name << " "; - } out << "<= " << rhs.toString(); } @@ -224,12 +230,12 @@ void SOSConstraint::violation(Rational& viol) const void SOSConstraint::print(std::ostream& out) const { out << name << " " << type; - if( sostype == TYPE_1 ) out << " 1: "; - else if( sostype == TYPE_2 ) out << " 2: "; + if( sostype == TYPE_1 ) + out << " 1: "; + else if( sostype == TYPE_2 ) + out << " 2: "; for( unsigned int i = 0; i < vars.size(); ++i ) - { out << vars[i]->name << " "; - } } bool SOSConstraint::checkType1(const Rational& tolerance) const @@ -242,7 +248,7 @@ bool SOSConstraint::checkType1(const Rational& tolerance) const // count number of non-zero variables for( unsigned int i = 0; i < vars.size(); ++i ) { - if((vars[i]->value < lb) || (vars[i]->value > ub)) + if( vars[i]->value < lb || vars[i]->value > ub ) cnt++; } return (cnt <= 1); @@ -260,10 +266,10 @@ bool SOSConstraint::checkType2(const Rational& tolerance) const // count number of non-zero variables for( unsigned int i = 0; i < vars.size(); ++i ) { - if((vars[i]->value < lb) || (vars[i]->value > ub)) + if( vars[i]->value < lb || vars[i]->value > ub ) { cnt++; - if( firstIndex == noIndex) + if( firstIndex == noIndex ) firstIndex = i; } } @@ -272,7 +278,7 @@ bool SOSConstraint::checkType2(const Rational& tolerance) const if( cnt < 2 ) return true; // check if var in position (firstIndex + 1) is non-zero - if((vars[firstIndex + 1]->value < lb) || (vars[firstIndex + 1]->value > ub)) + if( vars[firstIndex + 1]->value < lb || vars[firstIndex + 1]->value > ub ) return true; return false; } @@ -292,8 +298,10 @@ IndicatorConstraint::~IndicatorConstraint() bool IndicatorConstraint::check(const Rational& tolerance) const { Rational half(1,2); - if (ifvar->value > half && !ifvalue) return true; - if (ifvar->value < half && ifvalue) return true; + if( ifvar->value > half && !ifvalue ) + return true; + if( ifvar->value < half && ifvalue ) + return true; // thencons must be satisifed return thencons->check(tolerance); } @@ -302,8 +310,10 @@ void IndicatorConstraint::violation(Rational& viol) const { viol.toZero(); Rational half(1,2); - if (ifvar->value > half && !ifvalue) return; - if (ifvar->value < half && ifvalue) return; + if( ifvar->value > half && !ifvalue ) + return; + if( ifvar->value < half && ifvalue ) + return; // return violation of thencons thencons->violation(viol); } @@ -342,26 +352,28 @@ Model::~Model() Var* Model::getVar(const char* name) const { std::map::const_iterator itr = vars.find(name); - if( itr != vars.end() ) return itr->second; + if( itr != vars.end() ) + return itr->second; return NULL; } Constraint* Model::getCons(const char* name) const { std::map::const_iterator itr = conss.find(name); - if( itr != conss.end() ) return itr->second; + if( itr != conss.end() ) + return itr->second; return NULL; } void Model::pushVar(Var* var) { - assert ( var != NULL ); + assert( var != NULL ); vars[var->name] = var; } void Model::pushCons(Constraint* cons) { - assert ( cons != NULL ); + assert( cons != NULL ); conss[cons->name] = cons; } @@ -393,20 +405,19 @@ bool Model::readSol(const char* filename) } hasObjectiveValue = false; - bool hasVarValue = false; bool isSolFeas = true; - while(true) + while( true ) { // clear buffer content memset((void*)buf, 0, SOL_MAX_LINELEN); - if (NULL == fgets(buf, sizeof(buf), fp)) + if( fgets(buf, sizeof(buf), fp) == NULL ) break; // Normalize white spaces in line unsigned int len = strlen(buf); for( unsigned int i = 0; i < len; i++ ) - if( (buf[i] == '\t') || (buf[i] == '\n') || (buf[i] == '\r') ) + if( buf[i] == '\t' || buf[i] == '\n' || buf[i] == '\r' ) buf[i] = BLANK; // tokenize @@ -415,6 +426,7 @@ bool Model::readSol(const char* filename) if( varname == NULL ) continue; + // detect infeasible declaration if( strcmp(varname, "=infeas=") == 0 ) { isSolFeas = false; @@ -424,7 +436,28 @@ bool Model::readSol(const char* filename) const char* valuep = strtok_r(NULL, " ", &nexttok); assert( valuep != NULL ); + if( strcmp(varname, "no") == 0 && strcmp(valuep, "solution") == 0 ) + { + isSolFeas = false; + break; + } + + // skip solution status + if( strcmp(varname, "solution") == 0 && strcmp(valuep, "status:") == 0 ) + continue; + + // recognize objective value + bool isObjVal = false; + if( strcmp(varname, "=obj=") == 0 ) + isObjVal = true; + else if( strcmp(varname, "objective") == 0 && strcmp(valuep, "value:") == 0 ) + { + isObjVal = true; + valuep = strtok_r(NULL, " ", &nexttok); + } + + if( isObjVal ) { // read objective value hasObjectiveValue = true; @@ -434,15 +467,19 @@ bool Model::readSol(const char* filename) { // read variable value Var* var = getVar(varname); - if( var == NULL ) std::cerr << "unexpected variable<" << varname << "> in solution file" << std::endl; - assert( var != NULL ); + + if( var == NULL ) + { + std::cerr << "unexpected variable <" << varname << "> in solution file" << std::endl; + isSolFeas = false; + break; + } + Rational value; value.fromString(valuep); var->value = value; - hasVarValue = true; } } - isSolFeas = isSolFeas && hasVarValue; fclose(fp); fp = NULL; @@ -491,8 +528,10 @@ void Model::check( Rational prod; mult(prod, vitr->second->objCoef, vitr->second->value); - if (prod.isPositive()) objValPlus += prod; - else objvalMinus += prod; + if( prod.isPositive() ) + objValPlus += prod; + else + objvalMinus += prod; ++vitr; } @@ -513,7 +552,7 @@ void Model::check( Rational diff; sub(diff, objVal, objectiveValue); diff.abs(); - if (diff > objtol) + if( diff > objtol ) { correctObj = false; std::cerr << std::setprecision(16) << "Failed objective value check: " @@ -585,14 +624,17 @@ void Model::checkWrtExact( { feasibility = true; objective = true; - if( solverStatus == "unknown" ) return; + if( solverStatus == "unknown" ) + return; if( solverStatus == "stopped" ) { - if( exactStatus == "infeasible" ) std::cerr << "Warning: exactStatus says infeasible but we have a solution" << std::endl; + if( exactStatus == "infeasible" ) + std::cerr << "Warning: exactStatus says infeasible but we have a solution" << std::endl; return; } // we are in the case solverStatus == "solved" - if( exactStatus == "infeasible" ) feasibility = false; + if( exactStatus == "infeasible" ) + feasibility = false; // now we can check the objective value if( hasObjectiveValue ) { @@ -603,7 +645,8 @@ void Model::checkWrtExact( std::map::const_iterator vend = vars.end(); while( vitr != vend ) { - if( vitr->second->type == Var::CONTINUOUS ) isMIP = true; + if( vitr->second->type == Var::CONTINUOUS ) + isMIP = true; ++vitr; } Rational diff; @@ -615,7 +658,8 @@ void Model::checkWrtExact( << std::endl; // if error is greater than 0.1 fail in any case Rational maxAbsDiff(1, 10); - if( diff > maxAbsDiff ) objective = false; + if( diff > maxAbsDiff ) + objective = false; else { // try a relative error of 10^-7 (but only if this is not a pure integer problem!) @@ -625,8 +669,10 @@ void Model::checkWrtExact( magnitude.abs(); max(magnitude, one, magnitude); mult(relErr, magnitude, relErr); - if( isMIP && (diff < magnitude) ) objective = true; - else objective = false; + if( isMIP && diff < magnitude ) + objective = true; + else + objective = false; } } if( !objective ) diff --git a/check/solchecker/src/model.h b/check/solchecker/src/model.h index c05af8ac95..2e87388f64 100644 --- a/check/solchecker/src/model.h +++ b/check/solchecker/src/model.h @@ -14,6 +14,7 @@ #include #include + /** * @brief Class representing a problem variable. * Holds global variable information (such as bounds and objective coefficent)