/* $Id: IdiSolve.cpp 1878 2012-08-30 15:43:19Z forrest $ */ // 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 "CoinHelperFunctions.hpp" #include "Idiot.hpp" #define FIT #ifdef FIT #define HISTORY 8 #else #define HISTORY 7 #endif #define NSOLVE HISTORY-1 static void solveSmall(int nsolve, double **aIn, double **a, double * b) { int i, j; /* copy */ for (i = 0; i < nsolve; i++) { for (j = 0; j < nsolve; j++) { a[i][j] = aIn[i][j]; } } for (i = 0; i < nsolve; i++) { /* update using all previous */ double diagonal; int j; for (j = i; j < nsolve; j++) { int k; double value = a[i][j]; for (k = 0; k < i; k++) { value -= a[k][i] * a[k][j]; } a[i][j] = value; } diagonal = a[i][i]; if (diagonal < 1.0e-20) { diagonal = 0.0; } else { diagonal = 1.0 / sqrt(diagonal); } a[i][i] = diagonal; for (j = i + 1; j < nsolve; j++) { a[i][j] *= diagonal; } } /* update */ for (i = 0; i < nsolve; i++) { int j; double value = b[i]; for (j = 0; j < i; j++) { value -= b[j] * a[j][i]; } value *= a[i][i]; b[i] = value; } for (i = nsolve - 1; i >= 0; i--) { int j; double value = b[i]; for (j = i + 1; j < nsolve; j++) { value -= b[j] * a[i][j]; } value *= a[i][i]; b[i] = value; } } IdiotResult Idiot::objval(int nrows, int ncols, double * rowsol , double * colsol, double * pi, double * /*djs*/, const double * cost , const double * /*rowlower*/, const double * rowupper, const double * /*lower*/, const double * /*upper*/, const double * elemnt, const int * row, const CoinBigIndex * columnStart, const int * length, int extraBlock, int * rowExtra, double * solExtra, double * elemExtra, double * /*upperExtra*/, double * costExtra, double weight) { IdiotResult result; double objvalue = 0.0; double sum1 = 0.0, sum2 = 0.0; int i; for (i = 0; i < nrows; i++) { rowsol[i] = -rowupper[i]; } for (i = 0; i < ncols; i++) { CoinBigIndex j; double value = colsol[i]; if (value) { objvalue += value * cost[i]; if (elemnt) { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; rowsol[irow] += elemnt[j] * value; } } else { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; rowsol[irow] += value; } } } } /* adjust to make as feasible as possible */ /* no */ if (extraBlock) { for (i = 0; i < extraBlock; i++) { double element = elemExtra[i]; int irow = rowExtra[i]; objvalue += solExtra[i] * costExtra[i]; rowsol[irow] += solExtra[i] * element; } } for (i = 0; i < nrows; i++) { double value = rowsol[i]; sum1 += fabs(value); sum2 += value * value; pi[i] = -2.0 * weight * value; } result.infeas = sum1; result.objval = objvalue; result.weighted = objvalue + weight * sum2; result.sumSquared = sum2; return result; } IdiotResult Idiot::IdiSolve( int nrows, int ncols, double * COIN_RESTRICT rowsol , double * COIN_RESTRICT colsol, double * COIN_RESTRICT pi, double * COIN_RESTRICT djs, const double * COIN_RESTRICT origcost , double * COIN_RESTRICT rowlower, double * COIN_RESTRICT rowupper, const double * COIN_RESTRICT lower, const double * COIN_RESTRICT upper, const double * COIN_RESTRICT elemnt, const int * row, const CoinBigIndex * columnStart, const int * length, double * COIN_RESTRICT lambda, int maxIts, double mu, double drop, double maxmin, double offset, int strategy, double djTol, double djExit, double djFlag, CoinThreadRandom * randomNumberGenerator) { IdiotResult result; int i, j, k, iter; double value = 0.0, objvalue = 0.0, weightedObj = 0.0; double tolerance = 1.0e-8; double * history[HISTORY+1]; int ncolx; int nChange; int extraBlock = 0; int * rowExtra = NULL; double * COIN_RESTRICT solExtra = NULL; double * COIN_RESTRICT elemExtra = NULL; double * COIN_RESTRICT upperExtra = NULL; double * COIN_RESTRICT costExtra = NULL; double * COIN_RESTRICT useCostExtra = NULL; double * COIN_RESTRICT saveExtra = NULL; double * COIN_RESTRICT cost = NULL; double saveValue = 1.0e30; double saveOffset = offset; double useOffset = offset; /*#define NULLVECTOR*/ #ifndef NULLVECTOR int nsolve = NSOLVE; #else int nsolve = NSOLVE + 1; /* allow for null vector */ #endif int nflagged; double * COIN_RESTRICT thetaX; double * COIN_RESTRICT djX; double * COIN_RESTRICT bX; double * COIN_RESTRICT vX; double ** aX; double **aworkX; double ** allsum; double * COIN_RESTRICT saveSol = 0; const double * COIN_RESTRICT useCost = cost; double bestSol = 1.0e60; double weight = 0.5 / mu; char * statusSave = new char[2*ncols]; char * statusWork = statusSave + ncols; #define DJTEST 5 double djSave[DJTEST]; double largestDj = 0.0; double smallestDj = 1.0e60; double maxDj = 0.0; int doFull = 0; #define SAVEHISTORY 10 #define EVERY (2*SAVEHISTORY) #define AFTER SAVEHISTORY*(HISTORY+1) #define DROP 5 double after = AFTER; double obj[DROP]; double kbad = 0, kgood = 0; if (strategy & 128) after = 999999; /* no acceleration at all */ for (i = 0; i < DROP; i++) { obj[i] = 1.0e70; } //#define FOUR_GOES 2 #ifdef FOUR_GOES double * COIN_RESTRICT pi2 = new double [3*nrows]; double * COIN_RESTRICT rowsol2 = new double [3*nrows]; double * COIN_RESTRICT piX[4]; double * COIN_RESTRICT rowsolX[4]; int startsX[2][5]; int nChangeX[4]; double maxDjX[4]; double objvalueX[4]; int nflaggedX[4]; piX[0]=pi; piX[1]=pi2; piX[2]=pi2+nrows; piX[3]=piX[2]+nrows; rowsolX[0]=rowsol; rowsolX[1]=rowsol2; rowsolX[2]=rowsol2+nrows; rowsolX[3]=rowsolX[2]+nrows; #endif allsum = new double * [nsolve]; aX = new double * [nsolve]; aworkX = new double * [nsolve]; thetaX = new double[nsolve]; vX = new double[nsolve]; bX = new double[nsolve]; djX = new double[nsolve]; allsum[0] = pi; for (i = 0; i < nsolve; i++) { if (i) allsum[i] = new double[nrows]; aX[i] = new double[nsolve]; aworkX[i] = new double[nsolve]; } /* check = rows */ for (i = 0; i < nrows; i++) { if (rowupper[i] - rowlower[i] > tolerance) { extraBlock++; } } cost = new double[ncols]; memset(rowsol, 0, nrows * sizeof(double)); for (i = 0; i < ncols; i++) { CoinBigIndex j; double value = origcost[i] * maxmin; double value2 = colsol[i]; if (elemnt) { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; value += elemnt[j] * lambda[irow]; rowsol[irow] += elemnt[j] * value2; } } else { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; value += lambda[irow]; rowsol[irow] += value2; } } cost[i] = value; } if (extraBlock) { rowExtra = new int[extraBlock]; solExtra = new double[extraBlock]; elemExtra = new double[extraBlock]; upperExtra = new double[extraBlock]; costExtra = new double[extraBlock]; saveExtra = new double[extraBlock]; extraBlock = 0; int nbad = 0; for (i = 0; i < nrows; i++) { if (rowupper[i] - rowlower[i] > tolerance) { double smaller, difference; double value; saveExtra[extraBlock] = rowupper[i]; if (fabs(rowupper[i]) > fabs(rowlower[i])) { smaller = rowlower[i]; value = -1.0; } else { smaller = rowupper[i]; value = 1.0; } if (fabs(smaller) > 1.0e10) { if (!nbad) COIN_DETAIL_PRINT(printf("Can't handle rows where both bounds >1.0e10 %d %g\n", i, smaller)); nbad++; if (rowupper[i] < 0.0 || rowlower[i] > 0.0) abort(); if (fabs(rowupper[i]) > fabs(rowlower[i])) { rowlower[i] = -0.9e10; smaller = rowlower[i]; value = -1.0; } else { rowupper[i] = 0.9e10; saveExtra[extraBlock] = rowupper[i]; smaller = rowupper[i]; value = 1.0; } } difference = rowupper[i] - rowlower[i]; difference = CoinMin(difference, 1.0e31); rowupper[i] = smaller; elemExtra[extraBlock] = value; solExtra[extraBlock] = (rowupper[i] - rowsol[i]) / value; if (solExtra[extraBlock] < 0.0) solExtra[extraBlock] = 0.0; if (solExtra[extraBlock] > difference) solExtra[extraBlock] = difference; costExtra[extraBlock] = lambda[i] * value; upperExtra[extraBlock] = difference; rowsol[i] += value * solExtra[extraBlock]; rowExtra[extraBlock++] = i; } } if (nbad) COIN_DETAIL_PRINT(printf("%d bad values - results may be wrong\n", nbad)); } for (i = 0; i < nrows; i++) { offset += lambda[i] * rowsol[i]; } if ((strategy & 256) != 0) { /* save best solution */ saveSol = new double[ncols]; CoinMemcpyN(colsol, ncols, saveSol); if (extraBlock) { useCostExtra = new double[extraBlock]; memset(useCostExtra, 0, extraBlock * sizeof(double)); } useCost = origcost; useOffset = saveOffset; } else { useCostExtra = costExtra; useCost = cost; useOffset = offset; } ncolx = ncols + extraBlock; for (i = 0; i < HISTORY + 1; i++) { history[i] = new double[ncolx]; } for (i = 0; i < DJTEST; i++) { djSave[i] = 1.0e30; } for (i = 0; i < ncols; i++) { if (upper[i] - lower[i]) { statusSave[i] = 0; } else { statusSave[i] = 1; } } // for two pass method int start[2]; int stop[2]; int direction = -1; start[0] = 0; stop[0] = ncols; start[1] = 0; stop[1] = 0; iter = 0; for (; iter < maxIts; iter++) { double sum1 = 0.0, sum2 = 0.0; double lastObj = 1.0e70; int good = 0, doScale = 0; if (strategy & 16) { int ii = iter / EVERY + 1; ii = ii * EVERY; if (iter > ii - HISTORY * 2 && (iter & 1) == 0) { double * COIN_RESTRICT x = history[HISTORY-1]; for (i = HISTORY - 1; i > 0; i--) { history[i] = history[i-1]; } history[0] = x; CoinMemcpyN(colsol, ncols, history[0]); CoinMemcpyN(solExtra, extraBlock, history[0] + ncols); } } if ((iter % SAVEHISTORY) == 0 || doFull) { if ((strategy & 16) == 0) { double * COIN_RESTRICT x = history[HISTORY-1]; for (i = HISTORY - 1; i > 0; i--) { history[i] = history[i-1]; } history[0] = x; CoinMemcpyN(colsol, ncols, history[0]); CoinMemcpyN(solExtra, extraBlock, history[0] + ncols); } } /* start full try */ if ((iter % EVERY) == 0 || doFull) { // for next pass direction = - direction; // randomize. // The cast is to avoid gcc compiler warning int kcol = static_cast(ncols * randomNumberGenerator->randomDouble()); if (kcol == ncols) kcol = ncols - 1; if (direction > 0) { start[0] = kcol; stop[0] = ncols; start[1] = 0; stop[1] = kcol; #ifdef FOUR_GOES for (int itry=0;itry<2;itry++) { int chunk=(stop[itry]-start[itry]+FOUR_GOES-1)/FOUR_GOES; startsX[itry][0]=start[itry]; for (int i=1;i<5;i++) startsX[itry][i]=CoinMin(stop[itry],startsX[itry][i-1]+chunk); } #endif } else { start[0] = kcol; stop[0] = -1; start[1] = ncols - 1; stop[1] = kcol; #ifdef FOUR_GOES for (int itry=0;itry<2;itry++) { int chunk=(start[itry]-stop[itry]+FOUR_GOES-1)/FOUR_GOES; startsX[itry][0]=start[itry]; for (int i=1;i<5;i++) startsX[itry][i]=CoinMax(stop[itry],startsX[itry][i-1]-chunk); } #endif } int itry = 0; /*if ((strategy&16)==0) { double * COIN_RESTRICT x=history[HISTORY-1]; for (i=HISTORY-1;i>0;i--) { history[i]=history[i-1]; } history[0]=x; CoinMemcpyN(colsol,ncols,history[0]); CoinMemcpyN(solExtra,extraBlock,history[0]+ncols); }*/ while (!good) { itry++; #define MAXTRY 5 if (iter > after && doScale < 2 && itry < MAXTRY) { /* now full one */ for (i = 0; i < nrows; i++) { rowsol[i] = -rowupper[i]; } sum2 = 0.0; objvalue = 0.0; memset(pi, 0, nrows * sizeof(double)); { double * COIN_RESTRICT theta = thetaX; double * COIN_RESTRICT dj = djX; double * COIN_RESTRICT b = bX; double ** a = aX; double ** awork = aworkX; double * COIN_RESTRICT v = vX; double c; #ifdef FIT int ntot = 0, nsign = 0, ngood = 0, mgood[4] = {0, 0, 0, 0}; double diff1, diff2, val0, val1, val2, newValue; CoinMemcpyN(colsol, ncols, history[HISTORY-1]); CoinMemcpyN(solExtra, extraBlock, history[HISTORY-1] + ncols); #endif dj[0] = 0.0; for (i = 1; i < nsolve; i++) { dj[i] = 0.0; memset(allsum[i], 0, nrows * sizeof(double)); } for (i = 0; i < ncols; i++) { double value2 = colsol[i]; if (value2 > lower[i] + tolerance) { if(value2 < (upper[i] - tolerance)) { int k; objvalue += value2 * cost[i]; #ifdef FIT ntot++; val0 = history[0][i]; val1 = history[1][i]; val2 = history[2][i]; diff1 = val0 - val1; diff2 = val1 - val2; if (diff1*diff2 >= 0.0) { nsign++; if (fabs(diff1) < fabs(diff2)) { int ii = static_cast(fabs(4.0 * diff1 / diff2)); if (ii == 4) ii = 3; mgood[ii]++; ngood++; } if (fabs(diff1) < 0.75 * fabs(diff2)) { newValue = val1 + (diff1 * diff2) / (diff2 - diff1); } else { newValue = val1 + 4.0 * diff1; } } else { newValue = 0.333333333 * (val0 + val1 + val2); } if (newValue > upper[i] - tolerance) { newValue = upper[i]; } else if (newValue < lower[i] + tolerance) { newValue = lower[i]; } history[HISTORY-1][i] = newValue; #endif for (k = 0; k < HISTORY - 1; k++) { value = history[k][i] - history[k+1][i]; dj[k] += value * cost[i]; v[k] = value; } if (elemnt) { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; for (k = 0; k < HISTORY - 1; k++) { allsum[k][irow] += elemnt[j] * v[k]; } rowsol[irow] += elemnt[j] * value2; } } else { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; for (k = 0; k < HISTORY - 1; k++) { allsum[k][irow] += v[k]; } rowsol[irow] += value2; } } } else { /* at ub */ colsol[i] = upper[i]; value2 = colsol[i]; objvalue += value2 * cost[i]; if (elemnt) { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; rowsol[irow] += elemnt[j] * value2; } } else { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; rowsol[irow] += value2; } } } } else { /* at lb */ if (value2) { objvalue += value2 * cost[i]; if (elemnt) { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; rowsol[irow] += elemnt[j] * value2; } } else { for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; rowsol[irow] += value2; } } } } } #ifdef FIT /*printf("total %d, same sign %d, good %d %d %d %d %d\n", ntot,nsign,ngood,mgood[0],mgood[1],mgood[2],mgood[3]);*/ #endif if (extraBlock) { for (i = 0; i < extraBlock; i++) { double element = elemExtra[i]; int irow = rowExtra[i]; objvalue += solExtra[i] * costExtra[i]; if (solExtra[i] > tolerance && solExtra[i] < (upperExtra[i] - tolerance)) { double value2 = solExtra[i]; int k; for (k = 0; k < HISTORY - 1; k++) { value = history[k][i+ncols] - history[k+1][i+ncols]; dj[k] += value * costExtra[i]; allsum[k][irow] += element * value; } rowsol[irow] += element * value2; } else { double value2 = solExtra[i]; double element = elemExtra[i]; int irow = rowExtra[i]; rowsol[irow] += element * value2; } } } #ifdef NULLVECTOR if ((strategy & 64)) { double djVal = dj[0]; for (i = 0; i < ncols - nrows; i++) { double value2 = colsol[i]; if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) { value = history[0][i] - history[1][i]; } else { value = 0.0; } history[HISTORY][i] = value; } for (; i < ncols; i++) { double value2 = colsol[i]; double delta; int irow = i - (ncols - nrows); double oldSum = allsum[0][irow];; if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) { delta = history[0][i] - history[1][i]; } else { delta = 0.0; } djVal -= delta * cost[i]; oldSum -= delta; delta = - oldSum; djVal += delta * cost[i]; history[HISTORY][i] = delta; } dj[HISTORY-1] = djVal; djVal = 0.0; for (i = 0; i < ncols; i++) { double value2 = colsol[i]; if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance || i >= ncols - nrows) { int k; value = history[HISTORY][i]; djVal += value * cost[i]; for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) { int irow = row[j]; allsum[nsolve-1][irow] += value; } } } printf("djs %g %g\n", dj[HISTORY-1], djVal); } #endif for (i = 0; i < nsolve; i++) { int j; b[i] = 0.0; for (j = 0; j < nsolve; j++) { a[i][j] = 0.0; } } c = 0.0; for (i = 0; i < nrows; i++) { double value = rowsol[i]; for (k = 0; k < nsolve; k++) { v[k] = allsum[k][i]; b[k] += v[k] * value; } c += value * value; for (k = 0; k < nsolve; k++) { for (j = k; j < nsolve; j++) { a[k][j] += v[k] * v[j]; } } } sum2 = c; if (itry == 1) { lastObj = objvalue + weight * sum2; } for (k = 0; k < nsolve; k++) { b[k] = - (weight * b[k] + 0.5 * dj[k]); for (j = k; j < nsolve; j++) { a[k][j] *= weight; a[j][k] = a[k][j]; } } c *= weight; for (k = 0; k < nsolve; k++) { theta[k] = b[k]; } solveSmall(nsolve, a, awork, theta); if ((strategy & 64) != 0) { value = 10.0; for (k = 0; k < nsolve; k++) { value = CoinMax(value, fabs(theta[k])); } if (value > 10.0 && ((logLevel_ & 4) != 0)) { printf("theta %g %g %g\n", theta[0], theta[1], theta[2]); } value = 10.0 / value; for (k = 0; k < nsolve; k++) { theta[k] *= value; } } for (i = 0; i < ncolx; i++) { double valueh = 0.0; for (k = 0; k < HISTORY - 1; k++) { value = history[k][i] - history[k+1][i]; valueh += value * theta[k]; } #ifdef NULLVECTOR value = history[HISTORY][i]; valueh += value * theta[HISTORY-1]; #endif history[HISTORY][i] = valueh; } } #ifdef NULLVECTOR if ((strategy & 64)) { for (i = 0; i < ncols - nrows; i++) { if (colsol[i] <= lower[i] + tolerance || colsol[i] >= (upper[i] - tolerance)) { history[HISTORY][i] = 0.0;; } } tolerance = -tolerance; /* switch off test */ } #endif if (!doScale) { for (i = 0; i < ncols; i++) { if (colsol[i] > lower[i] + tolerance && colsol[i] < (upper[i] - tolerance)) { value = history[HISTORY][i]; colsol[i] += value; if (colsol[i] < lower[i] + tolerance) { colsol[i] = lower[i]; } else if (colsol[i] > upper[i] - tolerance) { colsol[i] = upper[i]; } } } if (extraBlock) { for (i = 0; i < extraBlock; i++) { if (solExtra[i] > tolerance && solExtra[i] < (upperExtra[i] - tolerance)) { value = history[HISTORY][i+ncols]; solExtra[i] += value; if (solExtra[i] < 0.0) { solExtra[i] = 0.0; } else if (solExtra[i] > upperExtra[i]) { solExtra[i] = upperExtra[i]; } } } } } else { double theta = 1.0; double saveTheta = theta; for (i = 0; i < ncols; i++) { if (colsol[i] > lower[i] + tolerance && colsol[i] < (upper[i] - tolerance)) { value = history[HISTORY][i]; if (value > 0) { if (theta * value + colsol[i] > upper[i]) { theta = (upper[i] - colsol[i]) / value; } } else if (value < 0) { if (colsol[i] + theta * value < lower[i]) { theta = (lower[i] - colsol[i]) / value; } } } } if (extraBlock) { for (i = 0; i < extraBlock; i++) { if (solExtra[i] > tolerance && solExtra[i] < (upperExtra[i] - tolerance)) { value = history[HISTORY][i+ncols]; if (value > 0) { if (theta * value + solExtra[i] > upperExtra[i]) { theta = (upperExtra[i] - solExtra[i]) / value; } } else if (value < 0) { if (solExtra[i] + theta * value < 0.0) { theta = -solExtra[i] / value; } } } } } if ((iter % 100 == 0) && (logLevel_ & 8) != 0) { if (theta < saveTheta) { printf(" - modified theta %g\n", theta); } } for (i = 0; i < ncols; i++) { if (colsol[i] > lower[i] + tolerance && colsol[i] < (upper[i] - tolerance)) { value = history[HISTORY][i]; colsol[i] += value * theta; } } if (extraBlock) { for (i = 0; i < extraBlock; i++) { if (solExtra[i] > tolerance && solExtra[i] < (upperExtra[i] - tolerance)) { value = history[HISTORY][i+ncols]; solExtra[i] += value * theta; } } } } #ifdef NULLVECTOR tolerance = fabs(tolerance); /* switch back on */ #endif if ((iter % 100) == 0 && (logLevel_ & 8) != 0) { printf("\n"); } } good = 1; result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost, rowlower, rowupper, lower, upper, elemnt, row, columnStart, length, extraBlock, rowExtra, solExtra, elemExtra, upperExtra, useCostExtra, weight); weightedObj = result.weighted; if (!iter) saveValue = weightedObj; objvalue = result.objval; sum1 = result.infeas; if (saveSol) { if (result.weighted < bestSol) { COIN_DETAIL_PRINT(printf("%d %g better than %g\n", iter, result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset)); bestSol = result.weighted; CoinMemcpyN(colsol, ncols, saveSol); } } #ifdef FITz if (iter > after) { IdiotResult result2; double ww, oo, ss; if (extraBlock) abort(); result2 = objval(nrows, ncols, row2, sol2, pi2, djs, cost, rowlower, rowupper, lower, upper, elemnt, row, columnStart, extraBlock, rowExtra, solExtra, elemExtra, upperExtra, costExtra, weight); ww = result2.weighted; oo = result2.objval; ss = result2.infeas; printf("wobj %g obj %g inf %g last %g\n", ww, oo, ss, lastObj); if (ww < weightedObj && ww < lastObj) { printf(" taken"); ntaken++; saving += weightedObj - ww; weightedObj = ww; objvalue = oo; sum1 = ss; CoinMemcpyN(row2, nrows, rowsol); CoinMemcpyN(pi2, nrows, pi); CoinMemcpyN(sol2, ncols, colsol); result = objval(nrows, ncols, rowsol, colsol, pi, djs, cost, rowlower, rowupper, lower, upper, elemnt, row, columnStart, extraBlock, rowExtra, solExtra, elemExtra, upperExtra, costExtra, weight); weightedObj = result.weighted; objvalue = result.objval; sum1 = result.infeas; if (ww < weightedObj) abort(); } else { printf(" not taken"); nottaken++; } } #endif /*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/ if (weightedObj > lastObj + 1.0e-4 && itry < MAXTRY) { if((logLevel_ & 16) != 0 && doScale) { printf("Weighted objective from %g to %g **** bad move\n", lastObj, weightedObj); } if (doScale) { good = 1; } if ((strategy & 3) == 1) { good = 0; if (weightedObj > lastObj + djExit) { if ((logLevel_ & 16) != 0) { printf("Weighted objective from %g to %g ?\n", lastObj, weightedObj); } CoinMemcpyN(history[0], ncols, colsol); CoinMemcpyN(history[0] + ncols, extraBlock, solExtra); good = 1; } } else if ((strategy & 3) == 2) { if (weightedObj > lastObj + 0.1 * maxDj) { CoinMemcpyN(history[0], ncols, colsol); CoinMemcpyN(history[0] + ncols, extraBlock, solExtra); doScale++; good = 0; } } else if ((strategy & 3) == 3) { if (weightedObj > lastObj + 0.001 * maxDj) { /*doScale++;*/ good = 0; } } } } if ((iter % checkFrequency_) == 0) { double best = weightedObj; double test = obj[0]; for (i = 1; i < DROP; i++) { obj[i-1] = obj[i]; if (best > obj[i]) best = obj[i]; } obj[DROP-1] = best; if (test - best < drop && (strategy & 8) == 0) { if ((logLevel_ & 8) != 0) { printf("Exiting as drop in %d its is %g after %d iterations\n", DROP * checkFrequency_, test - best, iter); } goto RETURN; } } if ((iter % logFreq_) == 0) { double piSum = 0.0; for (i = 0; i < nrows; i++) { piSum += (rowsol[i] + rowupper[i]) * pi[i]; } if ((logLevel_ & 2) != 0) { printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n", iter, sum1, objvalue * maxmin - useOffset, weightedObj - useOffset, piSum * maxmin - useOffset, maxDj); } } CoinMemcpyN(statusSave, ncols, statusWork); nflagged = 0; } nChange = 0; doFull = 0; maxDj = 0.0; // go through forwards or backwards and starting at odd places #ifdef FOUR_GOES for (int i=1;i>1; } else { istart=(start[itry]+stop[itry])>>1; istop = stop[itry]; } #endif #if 0 printf("istart %d istop %d direction %d array %d %d new %d %d\n", istart,istop,direction,start[itry],stop[itry], startsX[itry][iPar],startsX[itry][iPar+1]); #endif istart=startsX[itry][iPar]; istop=startsX[itry][iPar+1]; #else for (int itry = 0; itry < 2; itry++) { int istart = start[itry]; int istop = stop[itry]; #endif for (int icol=istart; icol != istop; icol += direction) { if (!statusWork[icol]) { CoinBigIndex j; double value = colsol[icol]; double djval = cost[icol]; double djval2, value2; double theta, a, b, c; if (elemnt) { for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) { int irow = row[j]; djval -= elemnt[j] * pi[irow]; } } else { for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) { int irow = row[j]; djval -= pi[irow]; } } /*printf("xx iter %d seq %d djval %g value %g\n", iter,i,djval,value);*/ if (djval > 1.0e-5) { value2 = (lower[icol] - value); } else { value2 = (upper[icol] - value); } djval2 = djval * value2; djval = fabs(djval); if (djval > djTol) { if (djval2 < -1.0e-4) { #ifndef FOUR_GOES nChange++; if (djval > maxDj) maxDj = djval; #else nChangeX[iPar]++; if (djval > maxDjX[iPar]) maxDjX[iPar] = djval; #endif /*if (djval>3.55e6) { printf("big\n"); }*/ a = 0.0; b = 0.0; c = 0.0; djval2 = cost[icol]; if (elemnt) { for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) { int irow = row[j]; double value = rowsol[irow]; c += value * value; a += elemnt[j] * elemnt[j]; b += value * elemnt[j]; } } else { for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) { int irow = row[j]; double value = rowsol[irow]; c += value * value; a += 1.0; b += value; } } a *= weight; b = b * weight + 0.5 * djval2; c *= weight; /* solve */ theta = -b / a; #ifndef FOUR_GOES if ((strategy & 4) != 0) { double valuep, thetap; value2 = a * theta * theta + 2.0 * b * theta; thetap = 2.0 * theta; valuep = a * thetap * thetap + 2.0 * b * thetap; if (valuep < value2 + djTol) { theta = thetap; kgood++; } else { kbad++; } } #endif if (theta > 0.0) { if (theta < upper[icol] - colsol[icol]) { value2 = theta; } else { value2 = upper[icol] - colsol[icol]; } } else { if (theta > lower[icol] - colsol[icol]) { value2 = theta; } else { value2 = lower[icol] - colsol[icol]; } } colsol[icol] += value2; #ifndef FOUR_GOES objvalue += cost[icol] * value2; #else objvalueX[iPar] += cost[icol] * value2; #endif if (elemnt) { for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) { int irow = row[j]; double value; rowsol[irow] += elemnt[j] * value2; value = rowsol[irow]; pi[irow] = -2.0 * weight * value; } } else { for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) { int irow = row[j]; double value; rowsol[irow] += value2; value = rowsol[irow]; pi[irow] = -2.0 * weight * value; } } } else { /* dj but at bound */ if (djval > djFlag) { statusWork[icol] = 1; #ifndef FOUR_GOES nflagged++; #else nflaggedX[iPar]++; #endif } } } } } #ifdef FOUR_GOES } #endif } #ifdef FOUR_GOES for (int i=0;i 1.0e-5) { value2 = -value; } else { value2 = (upperExtra[i] - value); } djval2 = djval * value2; if (djval2 < -1.0e-4 && fabs(djval) > djTol) { nChange++; a = 0.0; b = 0.0; c = 0.0; djval2 = costExtra[i]; value = rowsol[irow]; c += value * value; a += element * element; b += element * value; a *= weight; b = b * weight + 0.5 * djval2; c *= weight; /* solve */ theta = -b / a; if (theta > 0.0) { value2 = CoinMin(theta, upperExtra[i] - solExtra[i]); } else { value2 = CoinMax(theta, -solExtra[i]); } solExtra[i] += value2; rowsol[irow] += element * value2; value = rowsol[irow]; pi[irow] = -2.0 * weight * value; } } } if ((iter % 10) == 2) { for (int i = DJTEST - 1; i > 0; i--) { djSave[i] = djSave[i-1]; } djSave[0] = maxDj; largestDj = CoinMax(largestDj, maxDj); smallestDj = CoinMin(smallestDj, maxDj); for (int i = DJTEST - 1; i > 0; i--) { maxDj += djSave[i]; } maxDj = maxDj / static_cast (DJTEST); if (maxDj < djExit && iter > 50) { //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter); break; } if (nChange < 100) { djTol *= 0.5; } } } RETURN: if (kgood || kbad) { COIN_DETAIL_PRINT(printf("%g good %g bad\n", kgood, kbad)); } result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost, rowlower, rowupper, lower, upper, elemnt, row, columnStart, length, extraBlock, rowExtra, solExtra, elemExtra, upperExtra, useCostExtra, weight); result.djAtBeginning = largestDj; result.djAtEnd = smallestDj; result.dropThis = saveValue - result.weighted; if (saveSol) { if (result.weighted < bestSol) { bestSol = result.weighted; CoinMemcpyN(colsol, ncols, saveSol); } else { COIN_DETAIL_PRINT(printf("restoring previous - now %g best %g\n", result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset)); } } if (saveSol) { if (extraBlock) { delete [] useCostExtra; } CoinMemcpyN(saveSol, ncols, colsol); delete [] saveSol; } for (i = 0; i < nsolve; i++) { if (i) delete [] allsum[i]; delete [] aX[i]; delete [] aworkX[i]; } delete [] thetaX; delete [] djX; delete [] bX; delete [] vX; delete [] aX; delete [] aworkX; delete [] allsum; delete [] cost; #ifdef FOUR_GOES delete [] pi2 ; delete [] rowsol2 ; #endif for (i = 0; i < HISTORY + 1; i++) { delete [] history[i]; } delete [] statusSave; /* do original costs objvalue*/ result.objval = 0.0; for (i = 0; i < ncols; i++) { result.objval += colsol[i] * origcost[i]; } if (extraBlock) { for (i = 0; i < extraBlock; i++) { int irow = rowExtra[i]; rowupper[irow] = saveExtra[i]; } delete [] rowExtra; delete [] solExtra; delete [] elemExtra; delete [] upperExtra; delete [] costExtra; delete [] saveExtra; } result.iteration = iter; result.objval -= saveOffset; result.weighted = result.objval + weight * result.sumSquared; return result; }