/* $Id: Idiot.cpp 1931 2013-04-06 20:44:29Z stefan $ */ // Copyright (C) 2002, International Business Machines // Corporation and others. All Rights Reserved. // This code is licensed under the terms of the Eclipse Public License (EPL). #include "CoinPragma.hpp" #include #include #include #include #include "ClpPresolve.hpp" #include "Idiot.hpp" #include "CoinTime.hpp" #include "CoinSort.hpp" #include "CoinMessageHandler.hpp" #include "CoinHelperFunctions.hpp" #include "AbcCommon.hpp" // Redefine stuff for Clp #ifndef OSI_IDIOT #include "ClpMessage.hpp" #define OsiObjOffset ClpObjOffset #endif /**** strategy 4 - drop, exitDrop and djTolerance all relative: For first two major iterations these are small. Then: drop - exit a major iteration if drop over 5*checkFrequency < this is used as info->drop*(10.0+fabs(last weighted objective)) exitDrop - exit idiot if feasible and drop < this is used as info->exitDrop*(10.0+fabs(last objective)) djExit - exit a major iteration if largest dj (averaged over 5 checks) drops below this - used as info->djTolerance*(10.0+fabs(last weighted objective) djFlag - mostly skip variables with bad dj worse than this => 2*djExit djTol - only look at variables with dj better than this => 0.01*djExit ****************/ #define IDIOT_FIX_TOLERANCE 1e-6 #define SMALL_IDIOT_FIX_TOLERANCE 1e-10 int Idiot::dropping(IdiotResult result, double tolerance, double small, int *nbad) { if (result.infeas <= small) { double value = CoinMax(fabs(result.objval), fabs(result.dropThis)) + 1.0; if (result.dropThis > tolerance * value) { *nbad = 0; return 1; } else { (*nbad)++; if (*nbad > 4) { return 0; } else { return 1; } } } else { *nbad = 0; return 1; } } // Deals with whenUsed and slacks int Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd, double * colsol, const double * lower, const double * upper, const double * rowLower, const double * rowUpper, const double * cost, const double * element, double fixTolerance, double & objValue, double & infValue) { int n = 0; if ((strategy_ & 16384) == 0) { for (int i = ordinaryStart; i < ordinaryEnd; i++) { if (colsol[i] > lower[i] + fixTolerance) { if (colsol[i] < upper[i] - fixTolerance) { n++; } else { colsol[i] = upper[i]; } whenUsed_[i] = iteration; } else { colsol[i] = lower[i]; } } return n; } else { #ifdef COIN_DEVELOP printf("entering inf %g, obj %g\n", infValue, objValue); #endif int nrows = model_->getNumRows(); int ncols = model_->getNumCols(); int * posSlack = whenUsed_ + ncols; int * negSlack = posSlack + nrows; int * nextSlack = negSlack + nrows; double * rowsol = reinterpret_cast (nextSlack + ncols); memset(rowsol, 0, nrows * sizeof(double)); #ifdef OSI_IDIOT const CoinPackedMatrix * matrix = model_->getMatrixByCol(); #else ClpMatrixBase * matrix = model_->clpMatrix(); #endif const int * row = matrix->getIndices(); const CoinBigIndex * columnStart = matrix->getVectorStarts(); const int * columnLength = matrix->getVectorLengths(); //const double * element = matrix->getElements(); int i; objValue = 0.0; infValue = 0.0; for ( i = 0; i < ncols; i++) { if (nextSlack[i] == -1) { // not a slack if (colsol[i] > lower[i] + fixTolerance) { if (colsol[i] < upper[i] - fixTolerance) { n++; whenUsed_[i] = iteration; } else { colsol[i] = upper[i]; } whenUsed_[i] = iteration; } else { colsol[i] = lower[i]; } double value = colsol[i]; if (value) { objValue += cost[i] * value; CoinBigIndex j; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int iRow = row[j]; rowsol[iRow] += value * element[j]; } } } } // temp fix for infinite lbs - just limit to -1000 for (i = 0; i < nrows; i++) { double rowSave = rowsol[i]; int iCol; iCol = posSlack[i]; if (iCol >= 0) { // slide all slack down double rowValue = rowsol[i]; CoinBigIndex j = columnStart[iCol]; double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]); rowSave += (colsol[iCol] - lowerValue) * element[j]; colsol[iCol] = lowerValue; while (nextSlack[iCol] >= 0) { iCol = nextSlack[iCol]; double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]); j = columnStart[iCol]; rowSave += (colsol[iCol] - lowerValue) * element[j]; colsol[iCol] = lowerValue; } iCol = posSlack[i]; while (rowValue < rowLower[i] && iCol >= 0) { // want to increase double distance = rowLower[i] - rowValue; double value = element[columnStart[iCol]]; double thisCost = cost[iCol]; if (distance <= value*(upper[iCol] - colsol[iCol])) { // can get there double movement = distance / value; objValue += movement * thisCost; rowValue = rowLower[i]; colsol[iCol] += movement; } else { // can't get there double movement = upper[iCol] - colsol[iCol]; objValue += movement * thisCost; rowValue += movement * value; colsol[iCol] = upper[iCol]; iCol = nextSlack[iCol]; } } if (iCol >= 0) { // may want to carry on - because of cost? while (iCol >= 0 && cost[iCol] < 0 && rowValue < rowUpper[i]) { // want to increase double distance = rowUpper[i] - rowValue; double value = element[columnStart[iCol]]; double thisCost = cost[iCol]; if (distance <= value*(upper[iCol] - colsol[iCol])) { // can get there double movement = distance / value; objValue += movement * thisCost; rowValue = rowUpper[i]; colsol[iCol] += movement; iCol = -1; } else { // can't get there double movement = upper[iCol] - colsol[iCol]; objValue += movement * thisCost; rowValue += movement * value; colsol[iCol] = upper[iCol]; iCol = nextSlack[iCol]; } } if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance && colsol[iCol] < upper[iCol] - fixTolerance) { whenUsed_[i] = iteration; n++; } } rowsol[i] = rowValue; } iCol = negSlack[i]; if (iCol >= 0) { // slide all slack down double rowValue = rowsol[i]; CoinBigIndex j = columnStart[iCol]; double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]); rowSave += (colsol[iCol] - lowerValue) * element[j]; colsol[iCol] = lowerValue; while (nextSlack[iCol] >= 0) { iCol = nextSlack[iCol]; j = columnStart[iCol]; double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]); rowSave += (colsol[iCol] - lowerValue) * element[j]; colsol[iCol] = lowerValue; } iCol = negSlack[i]; while (rowValue > rowUpper[i] && iCol >= 0) { // want to increase double distance = -(rowUpper[i] - rowValue); double value = -element[columnStart[iCol]]; double thisCost = cost[iCol]; if (distance <= value*(upper[iCol] - lower[iCol])) { // can get there double movement = distance / value; objValue += movement * thisCost; rowValue = rowUpper[i]; colsol[iCol] += movement; } else { // can't get there double movement = upper[iCol] - lower[iCol]; objValue += movement * thisCost; rowValue -= movement * value; colsol[iCol] = upper[iCol]; iCol = nextSlack[iCol]; } } if (iCol >= 0) { // may want to carry on - because of cost? while (iCol >= 0 && cost[iCol] < 0 && rowValue > rowLower[i]) { // want to increase double distance = -(rowLower[i] - rowValue); double value = -element[columnStart[iCol]]; double thisCost = cost[iCol]; if (distance <= value*(upper[iCol] - colsol[iCol])) { // can get there double movement = distance / value; objValue += movement * thisCost; rowValue = rowLower[i]; colsol[iCol] += movement; iCol = -1; } else { // can't get there double movement = upper[iCol] - colsol[iCol]; objValue += movement * thisCost; rowValue -= movement * value; colsol[iCol] = upper[iCol]; iCol = nextSlack[iCol]; } } if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance && colsol[iCol] < upper[iCol] - fixTolerance) { whenUsed_[i] = iteration; n++; } } rowsol[i] = rowValue; } infValue += CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]); // just change rowsol[i] -= rowSave; } return n; } } /* returns -1 if none or start of costed slacks or -2 if there are costed slacks but it is messy */ static int countCostedSlacks(OsiSolverInterface * model) { #ifdef OSI_IDIOT const CoinPackedMatrix * matrix = model->getMatrixByCol(); #else ClpMatrixBase * matrix = model->clpMatrix(); #endif const int * row = matrix->getIndices(); const CoinBigIndex * columnStart = matrix->getVectorStarts(); const int * columnLength = matrix->getVectorLengths(); const double * element = matrix->getElements(); const double * rowupper = model->getRowUpper(); int nrows = model->getNumRows(); int ncols = model->getNumCols(); int slackStart = ncols - nrows; int nSlacks = nrows; int i; if (ncols <= nrows) return -1; while (1) { for (i = 0; i < nrows; i++) { int j = i + slackStart; CoinBigIndex k = columnStart[j]; if (columnLength[j] == 1) { if (row[k] != i || element[k] != 1.0) { nSlacks = 0; break; } } else { nSlacks = 0; break; } if (rowupper[i] <= 0.0) { nSlacks = 0; break; } } if (nSlacks || !slackStart) break; slackStart = 0; } if (!nSlacks) slackStart = -1; return slackStart; } void Idiot::crash(int numberPass, CoinMessageHandler * handler, const CoinMessages *messages, bool doCrossover) { // lightweight options int numberColumns = model_->getNumCols(); const double * objective = model_->getObjCoefficients(); int nnzero = 0; double sum = 0.0; int i; for (i = 0; i < numberColumns; i++) { if (objective[i]) { sum += fabs(objective[i]); nnzero++; } } sum /= static_cast (nnzero + 1); if (maxIts_ == 5) maxIts_ = 2; if (numberPass <= 0) majorIterations_ = static_cast(2 + log10(static_cast(numberColumns + 1))); else majorIterations_ = numberPass; // If mu not changed then compute if (mu_ == 1e-4) mu_ = CoinMax(1.0e-3, sum * 1.0e-5); if (maxIts2_ == 100) { if (!lightWeight_) { maxIts2_ = 105; } else if (lightWeight_ == 1) { mu_ *= 1000.0; maxIts2_ = 23; } else if (lightWeight_ == 2) { maxIts2_ = 11; } else { maxIts2_ = 23; } } //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_); solve2(handler, messages); #ifndef OSI_IDIOT if (doCrossover) { double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast (model_->numberRows()); if ((averageInfeas < 0.01 && (strategy_ & 512) != 0) || (strategy_ & 8192) != 0) crossOver(16 + 1); else crossOver(majorIterations_ < 1000000 ? 3 : 2); } #endif } void Idiot::solve() { CoinMessages dummy; solve2(NULL, &dummy); } void Idiot::solve2(CoinMessageHandler * handler, const CoinMessages * messages) { int strategy = 0; double d2; int i, n; int allOnes = 1; int iteration = 0; int iterationTotal = 0; int nTry = 0; /* number of tries at same weight */ double fixTolerance = IDIOT_FIX_TOLERANCE; int maxBigIts = maxBigIts_; int maxIts = maxIts_; int logLevel = logLevel_; int saveMajorIterations = majorIterations_; majorIterations_ = majorIterations_ % 1000000; if (handler) { if (handler->logLevel() > 0 && handler->logLevel() < 3) logLevel = 1; else if (!handler->logLevel()) logLevel = 0; else logLevel = 7; } double djExit = djTolerance_; double djFlag = 1.0 + 100.0 * djExit; double djTol = 0.00001; double mu = mu_; double drop = drop_; int maxIts2 = maxIts2_; double factor = muFactor_; double smallInfeas = smallInfeas_; double reasonableInfeas = reasonableInfeas_; double stopMu = stopMu_; double maxmin, offset; double lastWeighted = 1.0e50; double exitDrop = exitDrop_; double fakeSmall = smallInfeas; double firstInfeas; int badIts = 0; int slackStart, ordStart, ordEnd; int checkIteration = 0; int lambdaIteration = 0; int belowReasonable = 0; /* set if ever gone below reasonable infeas */ double bestWeighted = 1.0e60; double bestFeasible = 1.0e60; /* best solution while feasible */ IdiotResult result, lastResult; int saveStrategy = strategy_; const int strategies[] = {0, 2, 128}; double saveLambdaScale = 0.0; if ((saveStrategy & 128) != 0) { fixTolerance = SMALL_IDIOT_FIX_TOLERANCE; } #ifdef OSI_IDIOT const CoinPackedMatrix * matrix = model_->getMatrixByCol(); #else ClpMatrixBase * matrix = model_->clpMatrix(); #endif const int * row = matrix->getIndices(); const CoinBigIndex * columnStart = matrix->getVectorStarts(); const int * columnLength = matrix->getVectorLengths(); const double * element = matrix->getElements(); int nrows = model_->getNumRows(); int ncols = model_->getNumCols(); double * rowsol, * colsol; double * pi, * dj; #ifndef OSI_IDIOT double * cost = model_->objective(); double * lower = model_->columnLower(); double * upper = model_->columnUpper(); #else double * cost = new double [ncols]; CoinMemcpyN( model_->getObjCoefficients(), ncols, cost); const double * lower = model_->getColLower(); const double * upper = model_->getColUpper(); #endif const double *elemXX; double * saveSol; double * rowupper = new double[nrows]; // not const as modified CoinMemcpyN(model_->getRowUpper(), nrows, rowupper); double * rowlower = new double[nrows]; // not const as modified CoinMemcpyN(model_->getRowLower(), nrows, rowlower); CoinThreadRandom * randomNumberGenerator = model_->randomNumberGenerator(); int * whenUsed; double * lambda; saveSol = new double[ncols]; lambda = new double [nrows]; rowsol = new double[nrows]; colsol = new double [ncols]; CoinMemcpyN(model_->getColSolution(), ncols, colsol); pi = new double[nrows]; dj = new double[ncols]; delete [] whenUsed_; bool oddSlacks = false; // See if any costed slacks int numberSlacks = 0; for (i = 0; i < ncols; i++) { if (columnLength[i] == 1) numberSlacks++; } if (!numberSlacks) { whenUsed_ = new int[ncols]; } else { #ifdef COIN_DEVELOP printf("%d slacks\n", numberSlacks); #endif oddSlacks = true; int extra = static_cast (nrows * sizeof(double) / sizeof(int)); whenUsed_ = new int[2*ncols+2*nrows+extra]; int * posSlack = whenUsed_ + ncols; int * negSlack = posSlack + nrows; int * nextSlack = negSlack + nrows; for (i = 0; i < nrows; i++) { posSlack[i] = -1; negSlack[i] = -1; } for (i = 0; i < ncols; i++) nextSlack[i] = -1; for (i = 0; i < ncols; i++) { if (columnLength[i] == 1) { CoinBigIndex j = columnStart[i]; int iRow = row[j]; if (element[j] > 0.0) { if (posSlack[iRow] == -1) { posSlack[iRow] = i; } else { int iCol = posSlack[iRow]; while (nextSlack[iCol] >= 0) iCol = nextSlack[iCol]; nextSlack[iCol] = i; } } else { if (negSlack[iRow] == -1) { negSlack[iRow] = i; } else { int iCol = negSlack[iRow]; while (nextSlack[iCol] >= 0) iCol = nextSlack[iCol]; nextSlack[iCol] = i; } } } } // now sort for (i = 0; i < nrows; i++) { int iCol; iCol = posSlack[i]; if (iCol >= 0) { CoinBigIndex j = columnStart[iCol]; #ifndef NDEBUG int iRow = row[j]; #endif assert (element[j] > 0.0); assert (iRow == i); dj[0] = cost[iCol] / element[j]; whenUsed_[0] = iCol; int n = 1; while (nextSlack[iCol] >= 0) { iCol = nextSlack[iCol]; CoinBigIndex j = columnStart[iCol]; #ifndef NDEBUG int iRow = row[j]; #endif assert (element[j] > 0.0); assert (iRow == i); dj[n] = cost[iCol] / element[j]; whenUsed_[n++] = iCol; } for (j = 0; j < n; j++) { int jCol = whenUsed_[j]; nextSlack[jCol] = -2; } CoinSort_2(dj, dj + n, whenUsed_); // put back iCol = whenUsed_[0]; posSlack[i] = iCol; for (j = 1; j < n; j++) { int jCol = whenUsed_[j]; nextSlack[iCol] = jCol; iCol = jCol; } } iCol = negSlack[i]; if (iCol >= 0) { CoinBigIndex j = columnStart[iCol]; #ifndef NDEBUG int iRow = row[j]; #endif assert (element[j] < 0.0); assert (iRow == i); dj[0] = -cost[iCol] / element[j]; whenUsed_[0] = iCol; int n = 1; while (nextSlack[iCol] >= 0) { iCol = nextSlack[iCol]; CoinBigIndex j = columnStart[iCol]; #ifndef NDEBUG int iRow = row[j]; #endif assert (element[j] < 0.0); assert (iRow == i); dj[n] = -cost[iCol] / element[j]; whenUsed_[n++] = iCol; } for (j = 0; j < n; j++) { int jCol = whenUsed_[j]; nextSlack[jCol] = -2; } CoinSort_2(dj, dj + n, whenUsed_); // put back iCol = whenUsed_[0]; negSlack[i] = iCol; for (j = 1; j < n; j++) { int jCol = whenUsed_[j]; nextSlack[iCol] = jCol; iCol = jCol; } } } } whenUsed = whenUsed_; if (model_->getObjSense() == -1.0) { maxmin = -1.0; } else { maxmin = 1.0; } model_->getDblParam(OsiObjOffset, offset); if (!maxIts2) maxIts2 = maxIts; strategy = strategy_; strategy &= 3; memset(lambda, 0, nrows * sizeof(double)); slackStart = countCostedSlacks(model_); if (slackStart >= 0) { COIN_DETAIL_PRINT(printf("This model has costed slacks\n")); if (slackStart) { ordStart = 0; ordEnd = slackStart; } else { ordStart = nrows; ordEnd = ncols; } } else { ordStart = 0; ordEnd = ncols; } if (offset && logLevel > 2) { printf("** Objective offset is %g\n", offset); } /* compute reasonable solution cost */ for (i = 0; i < nrows; i++) { rowsol[i] = 1.0e31; } for (i = 0; i < ncols; i++) { CoinBigIndex j; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { if (element[j] != 1.0) { allOnes = 0; break; } } } if (allOnes) { elemXX = NULL; } else { elemXX = element; } // Do scaling if wanted bool scaled = false; #ifndef OSI_IDIOT if ((strategy_ & 32) != 0 && !allOnes) { if (model_->scalingFlag() > 0) scaled = model_->clpMatrix()->scale(model_) == 0; if (scaled) { const double * rowScale = model_->rowScale(); const double * columnScale = model_->columnScale(); double * oldLower = lower; double * oldUpper = upper; double * oldCost = cost; lower = new double[ncols]; upper = new double[ncols]; cost = new double[ncols]; CoinMemcpyN(oldLower, ncols, lower); CoinMemcpyN(oldUpper, ncols, upper); CoinMemcpyN(oldCost, ncols, cost); int icol, irow; for (icol = 0; icol < ncols; icol++) { double multiplier = 1.0 / columnScale[icol]; if (lower[icol] > -1.0e50) lower[icol] *= multiplier; if (upper[icol] < 1.0e50) upper[icol] *= multiplier; colsol[icol] *= multiplier; cost[icol] *= columnScale[icol]; } CoinMemcpyN(model_->rowLower(), nrows, rowlower); for (irow = 0; irow < nrows; irow++) { double multiplier = rowScale[irow]; if (rowlower[irow] > -1.0e50) rowlower[irow] *= multiplier; if (rowupper[irow] < 1.0e50) rowupper[irow] *= multiplier; rowsol[irow] *= multiplier; } int length = columnStart[ncols-1] + columnLength[ncols-1]; double * elemYY = new double[length]; for (i = 0; i < ncols; i++) { CoinBigIndex j; double scale = columnScale[i]; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int irow = row[j]; elemYY[j] = element[j] * scale * rowScale[irow]; } } elemXX = elemYY; } } #endif for (i = 0; i < ncols; i++) { CoinBigIndex j; double dd = columnLength[i]; dd = cost[i] / dd; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int irow = row[j]; if (dd < rowsol[irow]) { rowsol[irow] = dd; } } } d2 = 0.0; for (i = 0; i < nrows; i++) { d2 += rowsol[i]; } d2 *= 2.0; /* for luck */ d2 = d2 / static_cast (4 * nrows + 8000); d2 *= 0.5; /* halve with more flexible method */ if (d2 < 5.0) d2 = 5.0; if (djExit == 0.0) { djExit = d2; } if ((saveStrategy & 4) != 0) { /* go to relative tolerances - first small */ djExit = 1.0e-10; djFlag = 1.0e-5; drop = 1.0e-10; } memset(whenUsed, 0, ncols * sizeof(int)); strategy = strategies[strategy]; if ((saveStrategy & 8) != 0) strategy |= 64; /* don't allow large theta */ CoinMemcpyN(colsol, ncols, saveSol); lastResult = IdiSolve(nrows, ncols, rowsol , colsol, pi, dj, cost, rowlower, rowupper, lower, upper, elemXX, row, columnStart, columnLength, lambda, 0, mu, drop, maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator); // update whenUsed_ n = cleanIteration(iteration, ordStart, ordEnd, colsol, lower, upper, rowlower, rowupper, cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas); if ((strategy_ & 16384) != 0) { int * posSlack = whenUsed_ + ncols; int * negSlack = posSlack + nrows; int * nextSlack = negSlack + nrows; double * rowsol2 = reinterpret_cast (nextSlack + ncols); for (i = 0; i < nrows; i++) rowsol[i] += rowsol2[i]; } if ((logLevel_ & 1) != 0) { #ifndef OSI_IDIOT if (!handler) { #endif printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", iteration, lastResult.infeas, lastResult.objval, mu, lastResult.iteration, n); #ifndef OSI_IDIOT } else { handler->message(CLP_IDIOT_ITERATION, *messages) << iteration << lastResult.infeas << lastResult.objval << mu << lastResult.iteration << n << CoinMessageEol; } #endif } int numberBaseTrys = 0; // for first time int numberAway = -1; iterationTotal = lastResult.iteration; firstInfeas = lastResult.infeas; if ((strategy_ & 1024) != 0) reasonableInfeas = 0.5 * firstInfeas; if (lastResult.infeas < reasonableInfeas) lastResult.infeas = reasonableInfeas; double keepinfeas = 1.0e31; double lastInfeas = 1.0e31; double bestInfeas = 1.0e31; while ((mu > stopMu && lastResult.infeas > smallInfeas) || (lastResult.infeas <= smallInfeas && dropping(lastResult, exitDrop, smallInfeas, &badIts)) || checkIteration < 2 || lambdaIteration < lambdaIterations_) { if (lastResult.infeas <= exitFeasibility_) break; iteration++; checkIteration++; if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) { bestFeasible = lastResult.objval; } if (lastResult.infeas + mu * lastResult.objval < bestWeighted) { bestWeighted = lastResult.objval + mu * lastResult.objval; } if ((saveStrategy & 4096)) strategy |= 256; if ((saveStrategy & 4) != 0 && iteration > 2) { /* go to relative tolerances */ double weighted = 10.0 + fabs(lastWeighted); djExit = djTolerance_ * weighted; djFlag = 2.0 * djExit; drop = drop_ * weighted; djTol = 0.01 * djExit; } CoinMemcpyN(colsol, ncols, saveSol); result = IdiSolve(nrows, ncols, rowsol , colsol, pi, dj, cost, rowlower, rowupper, lower, upper, elemXX, row, columnStart, columnLength, lambda, maxIts, mu, drop, maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator); n = cleanIteration(iteration, ordStart, ordEnd, colsol, lower, upper, rowlower, rowupper, cost, elemXX, fixTolerance, result.objval, result.infeas); if ((strategy_ & 16384) != 0) { int * posSlack = whenUsed_ + ncols; int * negSlack = posSlack + nrows; int * nextSlack = negSlack + nrows; double * rowsol2 = reinterpret_cast (nextSlack + ncols); for (i = 0; i < nrows; i++) rowsol[i] += rowsol2[i]; } if ((logLevel_ & 1) != 0) { #ifndef OSI_IDIOT if (!handler) { #endif printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", iteration, result.infeas, result.objval, mu, result.iteration, n); #ifndef OSI_IDIOT } else { handler->message(CLP_IDIOT_ITERATION, *messages) << iteration << result.infeas << result.objval << mu << result.iteration << n << CoinMessageEol; } #endif } if (iteration > 50 && n == numberAway ) { if((result.infeas < 1.0e-4 && majorIterations_<200)||result.infeas<1.0e-8) { #ifdef CLP_INVESTIGATE printf("infeas small %g\n", result.infeas); #endif break; // not much happening } } if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) { if (lastInfeas != bestInfeas && CoinMin(result.infeas, lastInfeas) > 0.95 * bestInfeas) majorIterations_ = CoinMin(majorIterations_, iteration); // not getting feasible } lastInfeas = result.infeas; numberAway = n; keepinfeas = result.infeas; lastWeighted = result.weighted; iterationTotal += result.iteration; if (iteration == 1) { if ((strategy_ & 1024) != 0 && mu < 1.0e-10) result.infeas = firstInfeas * 0.8; if (majorIterations_ >= 50 || dropEnoughFeasibility_ <= 0.0) result.infeas *= 0.8; if (result.infeas > firstInfeas * 0.9 && result.infeas > reasonableInfeas) { iteration--; if (majorIterations_ < 50) mu *= 1.0e-1; else mu *= 0.7; bestFeasible = 1.0e31; bestWeighted = 1.0e60; numberBaseTrys++; if (mu < 1.0e-30 || (numberBaseTrys > 10 && lightWeight_)) { // back to all slack basis lightWeight_ = 2; break; } CoinMemcpyN(saveSol, ncols, colsol); } else { maxIts = maxIts2; checkIteration = 0; if ((strategy_ & 1024) != 0) mu *= 1.0e-1; } } else { } bestInfeas = CoinMin(bestInfeas, result.infeas); if (majorIterations_>100&&majorIterations_<200) { if (iteration==majorIterations_-100) { // redo double muX=mu*10.0; bestInfeas=1.0e3; mu=muX; nTry=0; } } if (iteration) { /* this code is in to force it to terminate sometime */ double changeMu = factor; if ((saveStrategy & 64) != 0) { keepinfeas = 0.0; /* switch off ranga's increase */ fakeSmall = smallInfeas; } else { fakeSmall = -1.0; } saveLambdaScale = 0.0; if (result.infeas > reasonableInfeas || (nTry + 1 == maxBigIts && result.infeas > fakeSmall)) { if (result.infeas > lastResult.infeas*(1.0 - dropEnoughFeasibility_) || nTry + 1 == maxBigIts || (result.infeas > lastResult.infeas * 0.9 && result.weighted > lastResult.weighted - dropEnoughWeighted_ * CoinMax(fabs(lastResult.weighted), fabs(result.weighted)))) { mu *= changeMu; if ((saveStrategy & 32) != 0 && result.infeas < reasonableInfeas && 0) { reasonableInfeas = CoinMax(smallInfeas, reasonableInfeas * sqrt(changeMu)); COIN_DETAIL_PRINT(printf("reasonable infeas now %g\n", reasonableInfeas)); } result.weighted = 1.0e60; nTry = 0; bestFeasible = 1.0e31; bestWeighted = 1.0e60; checkIteration = 0; lambdaIteration = 0; #define LAMBDA #ifdef LAMBDA if ((saveStrategy & 2048) == 0) { memset(lambda, 0, nrows * sizeof(double)); } #else memset(lambda, 0, nrows * sizeof(double)); #endif } else { nTry++; } } else if (lambdaIterations_ >= 0) { /* update lambda */ double scale = 1.0 / mu; int i, nnz = 0; saveLambdaScale = scale; lambdaIteration++; if ((saveStrategy & 4) == 0) drop = drop_ / 50.0; if (lambdaIteration > 4 && (((lambdaIteration % 10) == 0 && smallInfeas < keepinfeas) || ((lambdaIteration % 5) == 0 && 1.5 * smallInfeas < keepinfeas))) { //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas); smallInfeas *= 1.5; } if ((saveStrategy & 2048) == 0) { for (i = 0; i < nrows; i++) { if (lambda[i]) nnz++; lambda[i] += scale * rowsol[i]; } } else { nnz = 1; #ifdef LAMBDA for (i = 0; i < nrows; i++) { lambda[i] += scale * rowsol[i]; } #else for (i = 0; i < nrows; i++) { lambda[i] = scale * rowsol[i]; } for (i = 0; i < ncols; i++) { CoinBigIndex j; double value = cost[i] * maxmin; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int irow = row[j]; value += element[j] * lambda[irow]; } cost[i] = value * maxmin; } for (i = 0; i < nrows; i++) { offset += lambda[i] * rowupper[i]; lambda[i] = 0.0; } #ifdef DEBUG printf("offset %g\n", offset); #endif model_->setDblParam(OsiObjOffset, offset); #endif } nTry++; if (!nnz) { bestFeasible = 1.0e32; bestWeighted = 1.0e60; checkIteration = 0; result.weighted = 1.0e31; } #ifdef DEBUG double trueCost = 0.0; for (i = 0; i < ncols; i++) { int j; trueCost += cost[i] * colsol[i]; } printf("True objective %g\n", trueCost - offset); #endif } else { nTry++; } lastResult = result; if (result.infeas < reasonableInfeas && !belowReasonable) { belowReasonable = 1; bestFeasible = 1.0e32; bestWeighted = 1.0e60; checkIteration = 0; result.weighted = 1.0e31; } } if (iteration >= majorIterations_) { // If not feasible and crash then dive dive dive if (mu > 1.0e-12 && result.infeas > 1.0 && majorIterations_ < 40) { mu = 1.0e-30; majorIterations_ = iteration + 1; stopMu = 0.0; } else { if (logLevel > 2) printf("Exiting due to number of major iterations\n"); break; } } } majorIterations_ = saveMajorIterations; #ifndef OSI_IDIOT if (scaled) { // Scale solution and free arrays const double * rowScale = model_->rowScale(); const double * columnScale = model_->columnScale(); int icol, irow; for (icol = 0; icol < ncols; icol++) { colsol[icol] *= columnScale[icol]; saveSol[icol] *= columnScale[icol]; dj[icol] /= columnScale[icol]; } for (irow = 0; irow < nrows; irow++) { rowsol[irow] /= rowScale[irow]; pi[irow] *= rowScale[irow]; } // Don't know why getting Microsoft problems #if defined (_MSC_VER) delete [] ( double *) elemXX; #else delete [] elemXX; #endif model_->setRowScale(NULL); model_->setColumnScale(NULL); delete [] lower; delete [] upper; delete [] cost; lower = model_->columnLower(); upper = model_->columnUpper(); cost = model_->objective(); //rowlower = model_->rowLower(); } #endif #define TRYTHIS #ifdef TRYTHIS if ((saveStrategy & 2048) != 0) { double offset; model_->getDblParam(OsiObjOffset, offset); for (i = 0; i < ncols; i++) { CoinBigIndex j; double djval = cost[i] * maxmin; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int irow = row[j]; djval -= element[j] * lambda[irow]; } cost[i] = djval; } for (i = 0; i < nrows; i++) { offset += lambda[i] * rowupper[i]; } model_->setDblParam(OsiObjOffset, offset); } #endif if (saveLambdaScale) { /* back off last update */ for (i = 0; i < nrows; i++) { lambda[i] -= saveLambdaScale * rowsol[i]; } } muAtExit_ = mu; // For last iteration make as feasible as possible if (oddSlacks) strategy_ |= 16384; // not scaled n = cleanIteration(iteration, ordStart, ordEnd, colsol, lower, upper, model_->rowLower(), model_->rowUpper(), cost, element, fixTolerance, lastResult.objval, lastResult.infeas); #if 0 if ((logLevel & 1) == 0 || (strategy_ & 16384) != 0) { printf( "%d - mu %g, infeasibility %g, objective %g, %d interior\n", iteration, mu, lastResult.infeas, lastResult.objval, n); } #endif #ifndef OSI_IDIOT model_->setSumPrimalInfeasibilities(lastResult.infeas); #endif // Put back more feasible solution double saveInfeas[] = {0.0, 0.0}; for (int iSol = 0; iSol < 3; iSol++) { const double * solution = iSol ? colsol : saveSol; if (iSol == 2 && saveInfeas[0] < saveInfeas[1]) { // put back best solution CoinMemcpyN(saveSol, ncols, colsol); } double large = 0.0; int i; memset(rowsol, 0, nrows * sizeof(double)); for (i = 0; i < ncols; i++) { CoinBigIndex j; double value = solution[i]; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int irow = row[j]; rowsol[irow] += element[j] * value; } } for (i = 0; i < nrows; i++) { if (rowsol[i] > rowupper[i]) { double diff = rowsol[i] - rowupper[i]; if (diff > large) large = diff; } else if (rowsol[i] < rowlower[i]) { double diff = rowlower[i] - rowsol[i]; if (diff > large) large = diff; } } if (iSol < 2) saveInfeas[iSol] = large; if (logLevel > 2) printf("largest infeasibility is %g\n", large); } /* subtract out lambda */ for (i = 0; i < nrows; i++) { pi[i] -= lambda[i]; } for (i = 0; i < ncols; i++) { CoinBigIndex j; double djval = cost[i] * maxmin; for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int irow = row[j]; djval -= element[j] * pi[irow]; } dj[i] = djval; } if ((strategy_ & 1024) != 0) { double ratio = static_cast (ncols) / static_cast (nrows); COIN_DETAIL_PRINT(printf("col/row ratio %g infeas ratio %g\n", ratio, lastResult.infeas / firstInfeas)); if (lastResult.infeas > 0.01 * firstInfeas * ratio) { strategy_ &= (~1024); COIN_DETAIL_PRINT(printf(" - layer off\n")); } else { COIN_DETAIL_PRINT(printf(" - layer on\n")); } } delete [] saveSol; delete [] lambda; // save solution // duals not much use - but save anyway #ifndef OSI_IDIOT CoinMemcpyN(rowsol, nrows, model_->primalRowSolution()); CoinMemcpyN(colsol, ncols, model_->primalColumnSolution()); CoinMemcpyN(pi, nrows, model_->dualRowSolution()); CoinMemcpyN(dj, ncols, model_->dualColumnSolution()); #else model_->setColSolution(colsol); model_->setRowPrice(pi); delete [] cost; #endif delete [] rowsol; delete [] colsol; delete [] pi; delete [] dj; delete [] rowlower; delete [] rowupper; return ; } #ifndef OSI_IDIOT void Idiot::crossOver(int mode) { if (lightWeight_ == 2) { // total failure model_->allSlackBasis(); return; } double fixTolerance = IDIOT_FIX_TOLERANCE; #ifdef COIN_DEVELOP double startTime = CoinCpuTime(); #endif ClpSimplex * saveModel = NULL; ClpMatrixBase * matrix = model_->clpMatrix(); const int * row = matrix->getIndices(); const CoinBigIndex * columnStart = matrix->getVectorStarts(); const int * columnLength = matrix->getVectorLengths(); const double * element = matrix->getElements(); const double * rowupper = model_->getRowUpper(); int nrows = model_->getNumRows(); int ncols = model_->getNumCols(); double * rowsol, * colsol; // different for Osi double * lower = model_->columnLower(); double * upper = model_->columnUpper(); const double * rowlower = model_->getRowLower(); int * whenUsed = whenUsed_; rowsol = model_->primalRowSolution(); colsol = model_->primalColumnSolution();; double * cost = model_->objective(); int slackEnd, ordStart, ordEnd; int slackStart = countCostedSlacks(model_); int addAll = mode & 7; int presolve = 0; double djTolerance = djTolerance_; if (djTolerance > 0.0 && djTolerance < 1.0) djTolerance = 1.0; int iteration; int i, n = 0; double ratio = 1.0; double objValue = 0.0; if ((strategy_ & 128) != 0) { fixTolerance = SMALL_IDIOT_FIX_TOLERANCE; } if ((mode & 16) != 0 && addAll < 3) presolve = 1; double * saveUpper = NULL; double * saveLower = NULL; double * saveRowUpper = NULL; double * saveRowLower = NULL; bool allowInfeasible = ((strategy_ & 8192) != 0) || (majorIterations_ > 1000000); if (addAll < 3) { saveUpper = new double [ncols]; saveLower = new double [ncols]; CoinMemcpyN(upper, ncols, saveUpper); CoinMemcpyN(lower, ncols, saveLower); if (allowInfeasible) { saveRowUpper = new double [nrows]; saveRowLower = new double [nrows]; CoinMemcpyN(rowupper, nrows, saveRowUpper); CoinMemcpyN(rowlower, nrows, saveRowLower); double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast (model_->numberRows()); fixTolerance = CoinMax(fixTolerance, 1.0e-5 * averageInfeas); } } if (slackStart >= 0) { slackEnd = slackStart + nrows; if (slackStart) { ordStart = 0; ordEnd = slackStart; } else { ordStart = nrows; ordEnd = ncols; } } else { slackEnd = slackStart; ordStart = 0; ordEnd = ncols; } /* get correct rowsol (without known slacks) */ memset(rowsol, 0, nrows * sizeof(double)); for (i = ordStart; i < ordEnd; i++) { CoinBigIndex j; double value = colsol[i]; if (value < lower[i] + fixTolerance) { value = lower[i]; colsol[i] = value; } for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { int irow = row[j]; rowsol[irow] += value * element[j]; } } if (slackStart >= 0) { for (i = 0; i < nrows; i++) { if (ratio * rowsol[i] > rowlower[i] && rowsol[i] > 1.0e-8) { ratio = rowlower[i] / rowsol[i]; } } for (i = 0; i < nrows; i++) { rowsol[i] *= ratio; } for (i = ordStart; i < ordEnd; i++) { double value = colsol[i] * ratio; colsol[i] = value; objValue += value * cost[i]; } for (i = 0; i < nrows; i++) { double value = rowlower[i] - rowsol[i]; colsol[i+slackStart] = value; objValue += value * cost[i+slackStart]; } COIN_DETAIL_PRINT(printf("New objective after scaling %g\n", objValue)); } #if 0 maybe put back - but just get feasible ? // If not many fixed then just exit int numberFixed = 0; for (i = ordStart; i < ordEnd; i++) { if (colsol[i] < lower[i] + fixTolerance) numberFixed++; else if (colsol[i] > upper[i] - fixTolerance) numberFixed++; } if (numberFixed < ncols / 2) { addAll = 3; presolve = 0; } #endif #ifdef FEB_TRY int savePerturbation = model_->perturbation(); int saveOptions = model_->specialOptions(); model_->setSpecialOptions(saveOptions | 8192); if (savePerturbation_ == 50) model_->setPerturbation(56); #endif model_->createStatus(); /* addAll 0 - chosen,all used, all 1 - chosen, all 2 - all 3 - do not do anything - maybe basis */ for (i = ordStart; i < ordEnd; i++) { if (addAll < 2) { if (colsol[i] < lower[i] + fixTolerance) { upper[i] = lower[i]; colsol[i] = lower[i]; } else if (colsol[i] > upper[i] - fixTolerance) { lower[i] = upper[i]; colsol[i] = upper[i]; } } model_->setColumnStatus(i, ClpSimplex::superBasic); } if ((strategy_ & 16384) != 0) { // put in basis int * posSlack = whenUsed_ + ncols; int * negSlack = posSlack + nrows; int * nextSlack = negSlack + nrows; /* Laci - try both ways - to see what works - you can change second part as much as you want */ #ifndef LACI_TRY // was #if 1 // Array for sorting out slack values double * ratio = new double [ncols]; int * which = new int [ncols]; for (i = 0; i < nrows; i++) { if (posSlack[i] >= 0 || negSlack[i] >= 0) { int iCol; int nPlus = 0; int nMinus = 0; bool possible = true; // Get sum double sum = 0.0; iCol = posSlack[i]; while (iCol >= 0) { double value = element[columnStart[iCol]]; sum += value * colsol[iCol]; if (lower[iCol]) { possible = false; break; } else { nPlus++; } iCol = nextSlack[iCol]; } iCol = negSlack[i]; while (iCol >= 0) { double value = -element[columnStart[iCol]]; sum -= value * colsol[iCol]; if (lower[iCol]) { possible = false; break; } else { nMinus++; } iCol = nextSlack[iCol]; } //printf("%d plus, %d minus",nPlus,nMinus); //printf("\n"); if ((rowsol[i] - rowlower[i] < 1.0e-7 || rowupper[i] - rowsol[i] < 1.0e-7) && nPlus + nMinus < 2) possible = false; if (possible) { // Amount contributed by other varaibles sum = rowsol[i] - sum; double lo = rowlower[i]; if (lo > -1.0e20) lo -= sum; double up = rowupper[i]; if (up < 1.0e20) up -= sum; //printf("row bounds %g %g\n",lo,up); if (0) { double sum = 0.0; double x = 0.0; for (int k = 0; k < ncols; k++) { CoinBigIndex j; double value = colsol[k]; x += value * cost[k]; for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) { int irow = row[j]; if (irow == i) sum += element[j] * value; } } printf("Before sum %g <= %g <= %g cost %.18g\n", rowlower[i], sum, rowupper[i], x); } // set all to zero iCol = posSlack[i]; while (iCol >= 0) { colsol[iCol] = 0.0; iCol = nextSlack[iCol]; } iCol = negSlack[i]; while (iCol >= 0) { colsol[iCol] = 0.0; iCol = nextSlack[iCol]; } { int iCol; iCol = posSlack[i]; while (iCol >= 0) { //printf("col %d el %g sol %g bounds %g %g cost %g\n", // iCol,element[columnStart[iCol]], // colsol[iCol],lower[iCol],upper[iCol],cost[iCol]); iCol = nextSlack[iCol]; } iCol = negSlack[i]; while (iCol >= 0) { //printf("col %d el %g sol %g bounds %g %g cost %g\n", // iCol,element[columnStart[iCol]], // colsol[iCol],lower[iCol],upper[iCol],cost[iCol]); iCol = nextSlack[iCol]; } } //printf("now what?\n"); int n = 0; bool basic = false; if (lo > 0.0) { // Add in positive iCol = posSlack[i]; while (iCol >= 0) { double value = element[columnStart[iCol]]; ratio[n] = cost[iCol] / value; which[n++] = iCol; iCol = nextSlack[iCol]; } CoinSort_2(ratio, ratio + n, which); for (int i = 0; i < n; i++) { iCol = which[i]; double value = element[columnStart[iCol]]; if (lo >= upper[iCol]*value) { value *= upper[iCol]; sum += value; lo -= value; colsol[iCol] = upper[iCol]; } else { value = lo / value; sum += lo; lo = 0.0; colsol[iCol] = value; model_->setColumnStatus(iCol, ClpSimplex::basic); basic = true; } if (lo < 1.0e-7) break; } } else if (up < 0.0) { // Use lo so coding is more similar lo = -up; // Add in negative iCol = negSlack[i]; while (iCol >= 0) { double value = -element[columnStart[iCol]]; ratio[n] = cost[iCol] / value; which[n++] = iCol; iCol = nextSlack[iCol]; } CoinSort_2(ratio, ratio + n, which); for (int i = 0; i < n; i++) { iCol = which[i]; double value = -element[columnStart[iCol]]; if (lo >= upper[iCol]*value) { value *= upper[iCol]; sum += value; lo -= value; colsol[iCol] = upper[iCol]; } else { value = lo / value; sum += lo; lo = 0.0; colsol[iCol] = value; model_->setColumnStatus(iCol, ClpSimplex::basic); basic = true; } if (lo < 1.0e-7) break; } } if (0) { double sum2 = 0.0; double x = 0.0; for (int k = 0; k < ncols; k++) { CoinBigIndex j; double value = colsol[k]; x += value * cost[k]; for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) { int irow = row[j]; if (irow == i) sum2 += element[j] * value; } } printf("after sum %g <= %g <= %g cost %.18g (sum = %g)\n", rowlower[i], sum2, rowupper[i], x, sum); } rowsol[i] = sum; if (basic) { if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i])) model_->setRowStatus(i, ClpSimplex::atLowerBound); else model_->setRowStatus(i, ClpSimplex::atUpperBound); } } else { int n = 0; int iCol; iCol = posSlack[i]; while (iCol >= 0) { if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) { model_->setColumnStatus(iCol, ClpSimplex::basic); n++; } iCol = nextSlack[iCol]; } iCol = negSlack[i]; while (iCol >= 0) { if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) { model_->setColumnStatus(iCol, ClpSimplex::basic); n++; } iCol = nextSlack[iCol]; } if (n) { if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i])) model_->setRowStatus(i, ClpSimplex::atLowerBound); else model_->setRowStatus(i, ClpSimplex::atUpperBound); #ifdef CLP_INVESTIGATE if (n > 1) printf("%d basic on row %d!\n", n, i); #endif } } } } delete [] ratio; delete [] which; #else for (i = 0; i < nrows; i++) { int n = 0; int iCol; iCol = posSlack[i]; while (iCol >= 0) { if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) { model_->setColumnStatus(iCol, ClpSimplex::basic); n++; } iCol = nextSlack[iCol]; } iCol = negSlack[i]; while (iCol >= 0) { if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) { model_->setColumnStatus(iCol, ClpSimplex::basic); n++; } iCol = nextSlack[iCol]; } if (n) { if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i])) model_->setRowStatus(i, ClpSimplex::atLowerBound); else model_->setRowStatus(i, ClpSimplex::atUpperBound); #ifdef CLP_INVESTIGATE if (n > 1) printf("%d basic on row %d!\n", n, i); #endif } } #endif } double maxmin; if (model_->getObjSense() == -1.0) { maxmin = -1.0; } else { maxmin = 1.0; } bool justValuesPass = majorIterations_ > 1000000; if (slackStart >= 0) { for (i = 0; i < nrows; i++) { model_->setRowStatus(i, ClpSimplex::superBasic); } for (i = slackStart; i < slackEnd; i++) { model_->setColumnStatus(i, ClpSimplex::basic); } } else { /* still try and put singletons rather than artificials in basis */ int ninbas = 0; for (i = 0; i < nrows; i++) { model_->setRowStatus(i, ClpSimplex::basic); } for (i = 0; i < ncols; i++) { if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) { CoinBigIndex j = columnStart[i]; double value = element[j]; int irow = row[j]; double rlo = rowlower[irow]; double rup = rowupper[irow]; double clo = lower[i]; double cup = upper[i]; double csol = colsol[i]; /* adjust towards feasibility */ double move = 0.0; if (rowsol[irow] > rup) { move = (rup - rowsol[irow]) / value; if (value > 0.0) { /* reduce */ if (csol + move < clo) move = CoinMin(0.0, clo - csol); } else { /* increase */ if (csol + move > cup) move = CoinMax(0.0, cup - csol); } } else if (rowsol[irow] < rlo) { move = (rlo - rowsol[irow]) / value; if (value > 0.0) { /* increase */ if (csol + move > cup) move = CoinMax(0.0, cup - csol); } else { /* reduce */ if (csol + move < clo) move = CoinMin(0.0, clo - csol); } } else { /* move to improve objective */ if (cost[i]*maxmin > 0.0) { if (value > 0.0) { move = (rlo - rowsol[irow]) / value; /* reduce */ if (csol + move < clo) move = CoinMin(0.0, clo - csol); } else { move = (rup - rowsol[irow]) / value; /* increase */ if (csol + move > cup) move = CoinMax(0.0, cup - csol); } } else if (cost[i]*maxmin < 0.0) { if (value > 0.0) { move = (rup - rowsol[irow]) / value; /* increase */ if (csol + move > cup) move = CoinMax(0.0, cup - csol); } else { move = (rlo - rowsol[irow]) / value; /* reduce */ if (csol + move < clo) move = CoinMin(0.0, clo - csol); } } } rowsol[irow] += move * value; colsol[i] += move; /* put in basis if row was artificial */ if (rup - rlo < 1.0e-7 && model_->getRowStatus(irow) == ClpSimplex::basic) { model_->setRowStatus(irow, ClpSimplex::superBasic); model_->setColumnStatus(i, ClpSimplex::basic); ninbas++; } } } /*printf("%d in basis\n",ninbas);*/ } bool wantVector = false; if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) { // See if original wanted vector ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix()); wantVector = clpMatrixO->wantsSpecialColumnCopy(); } if (addAll < 3) { ClpPresolve pinfo; if (presolve) { if (allowInfeasible) { // fix up so will be feasible double * rhs = new double[nrows]; memset(rhs, 0, nrows * sizeof(double)); model_->clpMatrix()->times(1.0, colsol, rhs); double * rowupper = model_->rowUpper(); double * rowlower = model_->rowLower(); saveRowUpper = CoinCopyOfArray(rowupper, nrows); saveRowLower = CoinCopyOfArray(rowlower, nrows); double sum = 0.0; for (i = 0; i < nrows; i++) { if (rhs[i] > rowupper[i]) { sum += rhs[i] - rowupper[i]; rowupper[i] = rhs[i]; } if (rhs[i] < rowlower[i]) { sum += rowlower[i] - rhs[i]; rowlower[i] = rhs[i]; } } COIN_DETAIL_PRINT(printf("sum of infeasibilities %g\n", sum)); delete [] rhs; } saveModel = model_; pinfo.setPresolveActions(pinfo.presolveActions() | 16384); model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5); } if (model_) { if (!wantVector) { //#define TWO_GOES #ifdef ABC_INHERIT #ifndef TWO_GOES model_->dealWithAbc(1,justValuesPass ? 2 : 1); #else model_->dealWithAbc(1,1 + 11); #endif #else #ifndef TWO_GOES model_->primal(justValuesPass ? 2 : 1); #else model_->primal(1 + 11); #endif #endif } else { ClpMatrixBase * matrix = model_->clpMatrix(); ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); assert (clpMatrix); clpMatrix->makeSpecialColumnCopy(); #ifdef ABC_INHERIT model_->dealWithAbc(1,1); #else model_->primal(1); #endif clpMatrix->releaseSpecialColumnCopy(); } if (presolve) { model_->primal(); pinfo.postsolve(true); delete model_; model_ = saveModel; saveModel = NULL; } } else { // not feasible addAll = 1; presolve = 0; model_ = saveModel; saveModel = NULL; if (justValuesPass) #ifdef ABC_INHERIT model_->dealWithAbc(1,2); #else model_->primal(2); #endif } if (allowInfeasible) { CoinMemcpyN(saveRowUpper, nrows, model_->rowUpper()); CoinMemcpyN(saveRowLower, nrows, model_->rowLower()); delete [] saveRowUpper; delete [] saveRowLower; saveRowUpper = NULL; saveRowLower = NULL; } if (addAll < 2) { n = 0; if (!addAll ) { /* could do scans to get a good number */ iteration = 1; for (i = ordStart; i < ordEnd; i++) { if (whenUsed[i] >= iteration) { if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) { n++; upper[i] = saveUpper[i]; lower[i] = saveLower[i]; } } } } else { for (i = ordStart; i < ordEnd; i++) { if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) { n++; upper[i] = saveUpper[i]; lower[i] = saveLower[i]; } } delete [] saveUpper; delete [] saveLower; saveUpper = NULL; saveLower = NULL; } #ifdef COIN_DEVELOP printf("Time so far %g, %d now added from previous iterations\n", CoinCpuTime() - startTime, n); #endif if (justValuesPass) return; if (addAll) presolve = 0; if (presolve) { saveModel = model_; model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5); } else { presolve = 0; } if (!wantVector) { #ifdef ABC_INHERIT model_->dealWithAbc(1,1); #else model_->primal(1); #endif } else { ClpMatrixBase * matrix = model_->clpMatrix(); ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); assert (clpMatrix); clpMatrix->makeSpecialColumnCopy(); #ifdef ABC_INHERIT model_->dealWithAbc(1,1); #else model_->primal(1); #endif clpMatrix->releaseSpecialColumnCopy(); } if (presolve) { pinfo.postsolve(true); delete model_; model_ = saveModel; saveModel = NULL; } if (!addAll) { n = 0; for (i = ordStart; i < ordEnd; i++) { if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) { n++; upper[i] = saveUpper[i]; lower[i] = saveLower[i]; } } delete [] saveUpper; delete [] saveLower; saveUpper = NULL; saveLower = NULL; #ifdef COIN_DEVELOP printf("Time so far %g, %d now added from previous iterations\n", CoinCpuTime() - startTime, n); #endif } if (presolve) { saveModel = model_; model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5); } else { presolve = 0; } if (!wantVector) { #ifdef ABC_INHERIT model_->dealWithAbc(1,1); #else model_->primal(1); #endif } else { ClpMatrixBase * matrix = model_->clpMatrix(); ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); assert (clpMatrix); clpMatrix->makeSpecialColumnCopy(); #ifdef ABC_INHERIT model_->dealWithAbc(1,1); #else model_->primal(1); #endif clpMatrix->releaseSpecialColumnCopy(); } if (presolve) { pinfo.postsolve(true); delete model_; model_ = saveModel; saveModel = NULL; } } #ifdef COIN_DEVELOP printf("Total time in crossover %g\n", CoinCpuTime() - startTime); #endif delete [] saveUpper; delete [] saveLower; } #ifdef FEB_TRY model_->setSpecialOptions(saveOptions); model_->setPerturbation(savePerturbation); #endif return ; } #endif /*****************************************************************************/ // Default contructor Idiot::Idiot() { model_ = NULL; maxBigIts_ = 3; maxIts_ = 5; logLevel_ = 1; logFreq_ = 100; maxIts2_ = 100; djTolerance_ = 1e-1; mu_ = 1e-4; drop_ = 5.0; exitDrop_ = -1.0e20; muFactor_ = 0.3333; stopMu_ = 1e-12; smallInfeas_ = 1e-1; reasonableInfeas_ = 1e2; muAtExit_ = 1.0e31; strategy_ = 8; lambdaIterations_ = 0; checkFrequency_ = 100; whenUsed_ = NULL; majorIterations_ = 30; exitFeasibility_ = -1.0; dropEnoughFeasibility_ = 0.02; dropEnoughWeighted_ = 0.01; // adjust double nrows = 10000.0; int baseIts = static_cast (sqrt(static_cast(nrows))); baseIts = baseIts / 10; baseIts *= 10; maxIts2_ = 200 + baseIts + 5; maxIts2_ = 100; reasonableInfeas_ = static_cast (nrows) * 0.05; lightWeight_ = 0; } // Constructor from model Idiot::Idiot(OsiSolverInterface &model) { model_ = & model; maxBigIts_ = 3; maxIts_ = 5; logLevel_ = 1; logFreq_ = 100; maxIts2_ = 100; djTolerance_ = 1e-1; mu_ = 1e-4; drop_ = 5.0; exitDrop_ = -1.0e20; muFactor_ = 0.3333; stopMu_ = 1e-12; smallInfeas_ = 1e-1; reasonableInfeas_ = 1e2; muAtExit_ = 1.0e31; strategy_ = 8; lambdaIterations_ = 0; checkFrequency_ = 100; whenUsed_ = NULL; majorIterations_ = 30; exitFeasibility_ = -1.0; dropEnoughFeasibility_ = 0.02; dropEnoughWeighted_ = 0.01; // adjust double nrows; if (model_) nrows = model_->getNumRows(); else nrows = 10000.0; int baseIts = static_cast (sqrt(static_cast(nrows))); baseIts = baseIts / 10; baseIts *= 10; maxIts2_ = 200 + baseIts + 5; maxIts2_ = 100; reasonableInfeas_ = static_cast (nrows) * 0.05; lightWeight_ = 0; } // Copy constructor. Idiot::Idiot(const Idiot &rhs) { model_ = rhs.model_; if (model_ && rhs.whenUsed_) { int numberColumns = model_->getNumCols(); whenUsed_ = new int [numberColumns]; CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_); } else { whenUsed_ = NULL; } djTolerance_ = rhs.djTolerance_; mu_ = rhs.mu_; drop_ = rhs.drop_; muFactor_ = rhs.muFactor_; stopMu_ = rhs.stopMu_; smallInfeas_ = rhs.smallInfeas_; reasonableInfeas_ = rhs.reasonableInfeas_; exitDrop_ = rhs.exitDrop_; muAtExit_ = rhs.muAtExit_; exitFeasibility_ = rhs.exitFeasibility_; dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_; dropEnoughWeighted_ = rhs.dropEnoughWeighted_; maxBigIts_ = rhs.maxBigIts_; maxIts_ = rhs.maxIts_; majorIterations_ = rhs.majorIterations_; logLevel_ = rhs.logLevel_; logFreq_ = rhs.logFreq_; checkFrequency_ = rhs.checkFrequency_; lambdaIterations_ = rhs.lambdaIterations_; maxIts2_ = rhs.maxIts2_; strategy_ = rhs.strategy_; lightWeight_ = rhs.lightWeight_; } // Assignment operator. This copies the data Idiot & Idiot::operator=(const Idiot & rhs) { if (this != &rhs) { delete [] whenUsed_; model_ = rhs.model_; if (model_ && rhs.whenUsed_) { int numberColumns = model_->getNumCols(); whenUsed_ = new int [numberColumns]; CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_); } else { whenUsed_ = NULL; } djTolerance_ = rhs.djTolerance_; mu_ = rhs.mu_; drop_ = rhs.drop_; muFactor_ = rhs.muFactor_; stopMu_ = rhs.stopMu_; smallInfeas_ = rhs.smallInfeas_; reasonableInfeas_ = rhs.reasonableInfeas_; exitDrop_ = rhs.exitDrop_; muAtExit_ = rhs.muAtExit_; exitFeasibility_ = rhs.exitFeasibility_; dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_; dropEnoughWeighted_ = rhs.dropEnoughWeighted_; maxBigIts_ = rhs.maxBigIts_; maxIts_ = rhs.maxIts_; majorIterations_ = rhs.majorIterations_; logLevel_ = rhs.logLevel_; logFreq_ = rhs.logFreq_; checkFrequency_ = rhs.checkFrequency_; lambdaIterations_ = rhs.lambdaIterations_; maxIts2_ = rhs.maxIts2_; strategy_ = rhs.strategy_; lightWeight_ = rhs.lightWeight_; } return *this; } Idiot::~Idiot() { delete [] whenUsed_; }