From ba29b579833ef299f3503b4a14b512136fff5ec4 Mon Sep 17 00:00:00 2001 From: John Forrest Date: Wed, 17 Jul 2024 15:51:13 +0100 Subject: [PATCH] fix odd accuracy bug with some cut generators --- src/Osi/OsiCuts.cpp | 178 +++++++++++++++++++++++++++++++++++++++++++- src/Osi/OsiCuts.hpp | 6 ++ 2 files changed, 180 insertions(+), 4 deletions(-) diff --git a/src/Osi/OsiCuts.cpp b/src/Osi/OsiCuts.cpp index de1115719..399987a29 100644 --- a/src/Osi/OsiCuts.cpp +++ b/src/Osi/OsiCuts.cpp @@ -11,6 +11,7 @@ #include #include "OsiCuts.hpp" +#include "CoinHelperFunctions.hpp" //------------------------------------------------------------------- // Default Constructor @@ -278,7 +279,8 @@ OsiCuts::const_iterator OsiCuts::const_iterator::operator++() } return *this; } - +static bool scaleCutIntegral(double* cutElem, int* cutIndex, int cutNz, + double& cutRhs, double maxdelta); /* Insert a row cut unless it is a duplicate (CoinAbsFltEq) returns true if inserted */ bool OsiCuts::insertIfNotDuplicate(OsiRowCut &rc, CoinAbsFltEq treatAsSame) @@ -317,9 +319,9 @@ bool OsiCuts::insertIfNotDuplicate(OsiRowCut &rc, CoinAbsFltEq treatAsSame) } if (notDuplicate) { OsiRowCut *newCutPtr = new OsiRowCut(); + newCutPtr->setRow(vector); newCutPtr->setLb(newLb); newCutPtr->setUb(newUb); - newCutPtr->setRow(vector); newCutPtr->setGloballyValid(rc.globallyValid()); newCutPtr->setEffectiveness(rc.effectiveness()); rowCutPtrs_.push_back(newCutPtr); @@ -372,6 +374,174 @@ void OsiCuts::insertIfNotDuplicate(OsiRowCut &rc, CoinRelFltEq treatAsSame) rowCutPtrs_.push_back(newCutPtr); } } +#define USE__RATIONAL 60 +#include "CoinRational.hpp" +#define SMALL_VALUE1 1.0e-14 +static long computeGcd(long a, long b) { + // This is the standard Euclidean algorithm for gcd + long remainder = 1; + // Make sure a<=b (will always remain so) + if (a > b) { + // Swap a and b + long temp = a; + a = b; + b = temp; + } + // If zero then gcd is nonzero + if (!a) { + if (b) { + return b; + } + else { + return 0; + } + } + while (remainder) { + remainder = b % a; + b = a; + a = remainder; + } + return b; +} /* computeGcd */ +static bool scaleCutIntegral(double* cutElem, int* cutIndex, int cutNz, + double& cutRhs, double maxdelta) { + long gcd, lcm; + double maxscale = 1000; + long maxdnom = USE__RATIONAL; + //long numerator = 0, denominator = 0; + // Initialize gcd and lcm + CoinRational r = CoinRational(cutRhs, maxdelta, maxdnom); + if (r.getNumerator() != 0){ + gcd = labs(r.getNumerator()); + lcm = r.getDenominator(); + } + else{ + return false; + } + for (int i = 0; i < cutNz; ++i) { + CoinRational r = CoinRational(cutElem[i], maxdelta, maxdnom); + if (r.getNumerator() != 0){ + gcd = computeGcd(gcd, r.getNumerator()); + lcm *= r.getDenominator()/(computeGcd(lcm,r.getDenominator())); + } + else{ + return false; + } + } + double scale = ((double)lcm)/((double)gcd); + if (fabs(scale) > maxscale) { + return false; + } + scale = fabs(scale); + // Looks like we have a good scaling factor; scale and return; + for (int i = 0; i < cutNz; ++i) { + double value = cutElem[i]*scale; + cutElem[i] = floor(value+0.5); + assert (fabs(cutElem[i]-value)<1.0e-9); + } + { + double value = cutRhs*scale; + cutRhs = floor(value+0.5); + assert (fabs(cutRhs-value)<1.0e-9); + } + return true; +} /* scaleCutIntegral */ +// Returns value - floor but allowing for small errors +inline double above_integer(double value) { + double value2=floor(value); + double value3=floor(value+0.5); + if (fabs(value3-value)<1.0e-9*(fabs(value3)+1.0)) + return 0.0; + return value-value2; +} +/* Insert a row cut unless it is a duplicate - cut may get sorted. + Duplicate is defined as CoinAbsFltEq says same + returns true if inserted. + Also tries to "integerize" cut and checks for accuracy */ +bool OsiCuts::insertIfNotDuplicateAndClean(OsiRowCut &rc, + int typeCut , + CoinAbsFltEq treatAsSame) +{ + double newLb = rc.lb(); + double newUb = rc.ub(); + CoinPackedVector vector = rc.row(); + int numberElements = vector.getNumElements(); + int *newIndices = vector.getIndices(); + double *newElements = vector.getElements(); + CoinSort_2(newIndices, newIndices + numberElements, newElements); + bool notDuplicate = true; + int numberRowCuts = sizeRowCuts(); + for (int i = 0; i < numberRowCuts; i++) { + const OsiRowCut *cutPtr = rowCutPtr(i); + if (cutPtr->row().getNumElements() != numberElements) + continue; + if (!treatAsSame(cutPtr->lb(), newLb)) + continue; + if (!treatAsSame(cutPtr->ub(), newUb)) + continue; + const CoinPackedVector *thisVector = &(cutPtr->row()); + const int *indices = thisVector->getIndices(); + const double *elements = thisVector->getElements(); + int j; + for (j = 0; j < numberElements; j++) { + if (indices[j] != newIndices[j]) + break; + if (!treatAsSame(elements[j], newElements[j])) + break; + } + if (j == numberElements) { + notDuplicate = false; + break; + } + } + if (notDuplicate) { + OsiRowCut *newCutPtr = new OsiRowCut(); + // scale + if (newLb<-1.0e30 || newUb > 1.0e30) { + double maxdelta = 1.0e-12; //param.getEPS(); + double rhs = newLb<-1.0e30 ? newUb : newLb; + bool goodScale = scaleCutIntegral(newElements,newIndices,numberElements, + rhs,maxdelta); + if (false) { + // relax rhs a tiny bit + rhs += (newLb<-1.0e30) ? 1.0e-7 : -1.0e-7; + //if (numberCoefficients>=10||true) { + //rhs += 1.0e-7*sumCoefficients+1.0e-8*numberCoefficients; + //} + } + if (goodScale) { + if (newLb<-1.0e30) + newUb = rhs; + else + newLb = rhs; + } else { + // scale + double rhs = fabs((newLb<-1.0e30) ? newUb : newLb); + double largest = 1.0e-100; + double smallest = 1.0e100; + for (int i=0;ismallest*1.0e8||rhs>smallest*1.0e8) { + //printf ("badly scaled cut - rhs %g els %g -> %g - type %d\n",rhs,smallest,largest, + // typeCut); + //if (rhs>1.0e5 && smallest < 1.0e-4) + if (typeCut==61) // CglTwoMir + return false; + } + } + newCutPtr->setRow(numberElements,newIndices,newElements); + } else { + newCutPtr->setRow(vector); + } + newCutPtr->setLb(newLb); + newCutPtr->setUb(newUb); + newCutPtr->setGloballyValid(rc.globallyValid()); + newCutPtr->setEffectiveness(rc.effectiveness()); + rowCutPtrs_.push_back(newCutPtr); + } + return notDuplicate; +} -/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 -*/ diff --git a/src/Osi/OsiCuts.hpp b/src/Osi/OsiCuts.hpp index a10d98e19..783ae6231 100644 --- a/src/Osi/OsiCuts.hpp +++ b/src/Osi/OsiCuts.hpp @@ -148,6 +148,12 @@ class OSILIB_EXPORT OsiCuts { /** \brief Insert a row cut unless it is a duplicate - cut may get sorted. Duplicate is defined as CoinRelFltEq says same*/ void insertIfNotDuplicate(OsiRowCut &rc, CoinRelFltEq treatAsSame); + /** \brief Insert a row cut unless it is a duplicate - cut may get sorted. + Duplicate is defined as CoinAbsFltEq says same + returns true if inserted. + Also tries to "integerize" cut and checks for accuracy */ + bool insertIfNotDuplicateAndClean(OsiRowCut &rc, + int cutType, CoinAbsFltEq treatAsSame = CoinAbsFltEq(1.0e-12)); /** \brief Insert a column cut */ inline void insert(const OsiColCut &cc);