Changeset 91 for trunk/CbcHeuristic.cpp


Ignore:
Timestamp:
Mar 28, 2005 3:53:37 PM (14 years ago)
Author:
forrest
Message:

allow more than singletons

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/CbcHeuristic.cpp

    r75 r91  
    216216  double penalty=0.0;
    217217 
    218   // see if feasible
     218  // see if feasible - just using singletons
    219219  for (i=0;i<numberRows;i++) {
    220220    double value = rowActivity[i];
     
    236236      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
    237237        int iColumn = column[k];
     238        // See if all elements help
    238239        if (columnLength[iColumn]==1) {
    239240          double currentValue = newSolution[iColumn];
     
    302303      }
    303304      penalty += fabs(thisInfeasibility);
     305    }
     306  }
     307  if (penalty) {
     308    // see if feasible using any
     309    // first continuous
     310    double penaltyChange=0.0;
     311    int iColumn;
     312    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     313      if (solver->isInteger(iColumn))
     314        continue;
     315      double currentValue = newSolution[iColumn];
     316      double lowerValue = lower[iColumn];
     317      double upperValue = upper[iColumn];
     318      int j;
     319      int anyBadDown=0;
     320      int anyBadUp=0;
     321      double upImprovement=0.0;
     322      double downImprovement=0.0;
     323      for (j=columnStart[iColumn];
     324           j<columnStart[iColumn]+columnLength[iColumn];j++) {
     325        int iRow = row[j];
     326        if (rowUpper[iRow]>rowLower[iRow]) {
     327          double value = element[j];
     328          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
     329            // infeasible above
     330            downImprovement += value;
     331            upImprovement -= value;
     332            if (value>0.0)
     333              anyBadUp++;
     334            else
     335              anyBadDown++;
     336          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
     337            // feasible at ub
     338            if (value>0.0) {
     339              upImprovement -= value;
     340              anyBadUp++;
     341            } else {
     342              downImprovement += value;
     343              anyBadDown++;
     344            }
     345          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
     346            // feasible in interior
     347          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
     348            // feasible at lb
     349            if (value<0.0) {
     350              upImprovement += value;
     351              anyBadUp++;
     352            } else {
     353              downImprovement -= value;
     354              anyBadDown++;
     355            }
     356          } else {
     357            // infeasible below
     358            downImprovement -= value;
     359            upImprovement += value;
     360            if (value<0.0)
     361              anyBadUp++;
     362            else
     363              anyBadDown++;
     364          }
     365        } else {
     366          // equality row
     367          double value = element[j];
     368          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
     369            // infeasible above
     370            downImprovement += value;
     371            upImprovement -= value;
     372            if (value>0.0)
     373              anyBadUp++;
     374            else
     375              anyBadDown++;
     376          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
     377            // infeasible below
     378            downImprovement -= value;
     379            upImprovement += value;
     380            if (value<0.0)
     381              anyBadUp++;
     382            else
     383              anyBadDown++;
     384          } else {
     385            // feasible - no good
     386            anyBadUp=-1;
     387            break;
     388          }
     389        }
     390      }
     391      // could change tests for anyBad
     392      if (anyBadUp)
     393        upImprovement=0.0;
     394      if (anyBadDown)
     395        downImprovement=0.0;
     396      double way=0.0;
     397      double improvement=0.0;
     398      if (downImprovement>0.0&&currentValue>lowerValue) {
     399        way=-1.0;
     400        improvement = downImprovement;
     401      } else if (upImprovement>0.0&&currentValue<upperValue) {
     402        way=1.0;
     403        improvement = upImprovement;
     404      }
     405      if (way) {
     406        // can improve
     407        double distance=COIN_DBL_MAX;
     408        for (j=columnStart[iColumn];
     409             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     410          int iRow = row[j];
     411          double value = element[j]*way;
     412          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
     413            // infeasible above
     414            assert (value<0.0);
     415            double gap = rowActivity[iRow]-rowUpper[iRow];
     416            if (gap+value*distance<0.0)
     417              distance = -gap/value;
     418          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
     419            // infeasible below
     420            assert (value>0.0);
     421            double gap = rowActivity[iRow]-rowLower[iRow];
     422            if (gap+value*distance>0.0)
     423              distance = -gap/value;
     424          } else {
     425            // feasible
     426            if (value>0) {
     427              double gap = rowActivity[iRow]-rowUpper[iRow];
     428              if (gap+value*distance>0.0)
     429              distance = -gap/value;
     430            } else {
     431              double gap = rowActivity[iRow]-rowLower[iRow];
     432              if (gap+value*distance<0.0)
     433                distance = -gap/value;
     434            }
     435          }
     436        }
     437        //move
     438        penaltyChange += improvement*distance;
     439        distance *= way;
     440        newSolution[iColumn] += distance;
     441        newSolutionValue += direction*objective[iColumn]*distance;
     442        for (j=columnStart[iColumn];
     443             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     444          int iRow = row[j];
     445          double value = element[j];
     446          rowActivity[iRow] += distance*value;
     447        }
     448      }
     449    }
     450    // and now all if improving
     451    double lastChange= penaltyChange ? 1.0 : 0.0;
     452    while (lastChange>1.0e-2) {
     453      lastChange=0;
     454      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     455        bool isInteger = solver->isInteger(iColumn);
     456        double currentValue = newSolution[iColumn];
     457        double lowerValue = lower[iColumn];
     458        double upperValue = upper[iColumn];
     459        int j;
     460        int anyBadDown=0;
     461        int anyBadUp=0;
     462        double upImprovement=0.0;
     463        double downImprovement=0.0;
     464        for (j=columnStart[iColumn];
     465             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     466          int iRow = row[j];
     467          double value = element[j];
     468          if (isInteger) {
     469            if (value>0.0) {
     470              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
     471                anyBadUp++;
     472              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
     473                anyBadDown++;
     474            } else {
     475              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
     476                anyBadDown++;
     477              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
     478                anyBadUp++;
     479            }
     480          }
     481          if (rowUpper[iRow]>rowLower[iRow]) {
     482            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
     483              // infeasible above
     484              downImprovement += value;
     485              upImprovement -= value;
     486              if (value>0.0)
     487                anyBadUp++;
     488              else
     489                anyBadDown++;
     490            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
     491              // feasible at ub
     492              if (value>0.0) {
     493                upImprovement -= value;
     494                anyBadUp++;
     495              } else {
     496                downImprovement += value;
     497                anyBadDown++;
     498              }
     499            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
     500              // feasible in interior
     501            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
     502              // feasible at lb
     503              if (value<0.0) {
     504                upImprovement += value;
     505                anyBadUp++;
     506              } else {
     507                downImprovement -= value;
     508                anyBadDown++;
     509              }
     510            } else {
     511              // infeasible below
     512              downImprovement -= value;
     513              upImprovement += value;
     514              if (value<0.0)
     515                anyBadUp++;
     516              else
     517                anyBadDown++;
     518            }
     519          } else {
     520            // equality row
     521            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
     522              // infeasible above
     523              downImprovement += value;
     524              upImprovement -= value;
     525              if (value>0.0)
     526                anyBadUp++;
     527              else
     528                anyBadDown++;
     529            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
     530              // infeasible below
     531              downImprovement -= value;
     532              upImprovement += value;
     533              if (value<0.0)
     534                anyBadUp++;
     535              else
     536                anyBadDown++;
     537            } else {
     538              // feasible - no good
     539              anyBadUp=-1;
     540              anyBadDown=-1;
     541              break;
     542            }
     543          }
     544        }
     545        // could change tests for anyBad
     546        if (anyBadUp)
     547          upImprovement=0.0;
     548        if (anyBadDown)
     549          downImprovement=0.0;
     550        double way=0.0;
     551        double improvement=0.0;
     552        if (downImprovement>0.0&&currentValue>lowerValue) {
     553          way=-1.0;
     554          improvement = downImprovement;
     555        } else if (upImprovement>0.0&&currentValue<upperValue) {
     556          way=1.0;
     557          improvement = upImprovement;
     558        }
     559        if (way) {
     560          // can improve
     561          double distance=COIN_DBL_MAX;
     562          for (j=columnStart[iColumn];
     563               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     564            int iRow = row[j];
     565            double value = element[j]*way;
     566            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
     567              // infeasible above
     568              assert (value<0.0);
     569              double gap = rowActivity[iRow]-rowUpper[iRow];
     570              if (gap+value*distance<0.0)
     571                distance = -gap/value;
     572            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
     573              // infeasible below
     574              assert (value>0.0);
     575              double gap = rowActivity[iRow]-rowLower[iRow];
     576              if (gap+value*distance>0.0)
     577                distance = -gap/value;
     578            } else {
     579              // feasible
     580              if (value>0) {
     581                double gap = rowActivity[iRow]-rowUpper[iRow];
     582                if (gap+value*distance>0.0)
     583                  distance = -gap/value;
     584              } else {
     585                double gap = rowActivity[iRow]-rowLower[iRow];
     586                if (gap+value*distance<0.0)
     587                  distance = -gap/value;
     588              }
     589            }
     590          }
     591          if (isInteger)
     592            distance = floor(distance+1.0e-8);
     593          if (!distance) {
     594            // should never happen
     595            printf("zero distance in CbcRounding - debug\n");
     596          }
     597          //move
     598          lastChange += improvement*distance;
     599          distance *= way;
     600          newSolution[iColumn] += distance;
     601          newSolutionValue += direction*objective[iColumn]*distance;
     602          for (j=columnStart[iColumn];
     603               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     604            int iRow = row[j];
     605            double value = element[j];
     606            rowActivity[iRow] += distance*value;
     607          }
     608        }
     609      }
     610      penaltyChange += lastChange;
     611    }
     612    penalty -= penaltyChange;
     613    if (penalty<1.0e-5*fabs(penaltyChange)) {
     614      // recompute
     615      penalty=0.0;
     616      for (i=0;i<numberRows;i++) {
     617        double value = rowActivity[i];
     618        if (value<rowLower[i]-primalTolerance)
     619          penalty += rowLower[i]-value;
     620        else if (value>rowUpper[i]+primalTolerance)
     621          penalty += value-rowUpper[i];
     622      }
    304623    }
    305624  }
Note: See TracChangeset for help on using the changeset viewer.