Changeset 191 for branches


Ignore:
Timestamp:
Aug 7, 2003 2:49:43 AM (16 years ago)
Author:
forrest
Message:

more stuff

Location:
branches/pre
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r190 r191  
    565565      // may factorize, checks if problem finished
    566566      statusOfProblemInPrimal(lastCleaned,factorType,progress_,info);
     567      if (problemStatus_>=0)
     568        break; // declare victory
    567569     
    568570      // Compute objective function from scratch
     
    11381140        }
    11391141        {
     1142          double oldValue = objectiveValue_;
    11401143          // Compute objective function from scratch
    11411144          const CoinPackedMatrix * quadratic =
     
    11921195          printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
    11931196                 objectiveValue_,infeasCost,objectiveValue_+infeasCost);
     1197          if (objectiveValue_>oldValue) {
     1198            // make less likely this one will come in again
     1199            double multiplier = CoinDrand48();
     1200            int iSequence;
     1201            if (!cleanupIteration) {
     1202              iSequence = sequenceIn_;
     1203            } else {
     1204              int iSequence = info->crucialSj();
     1205              if (iSequence>numberRows_)
     1206                iSequence -= numberRows_;
     1207              else
     1208                iSequence = numberColumns_ +(iSequence-numberXColumns);
     1209            }
     1210            info->djWeight()[iSequence] /= (5.0+multiplier);
     1211          }
    11941212        }
    11951213        if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
     
    12141232        cleanupIteration = true;
    12151233        // may not be correct on second time
    1216         if (result&&(returnCode==-1||returnCode==-2)) {
     1234        if (result&&(returnCode==-1||returnCode==-2||returnCode==-3)) {
    12171235          assert(sequenceOut_<numberXColumns||
    12181236                 sequenceOut_>=numberColumns_);
     
    15561574  double upperTheta = maximumMovement;
    15571575  bool throwAway=false;
    1558   if (numberIterations_==700||numberIterations_==694) {
     1576  if (numberIterations_>=2007) {
    15591577    printf("Bad iteration coming up after iteration %d\n",numberIterations_);
    15601578  }
     
    16451663      maximumSwapped = max(maximumSwapped,numberSwapped);
    16461664
    1647       double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 0.1*infeasibilityCost_;
    1648       // but make a bit more pessimistic
    1649       dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
     1665      //double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 0.1*infeasibilityCost_;
     1666      double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1;
     1667      // but make a bit more pessimistic if pivot
     1668      //if (bestPivot>acceptablePivot)
     1669      //dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
    16501670      //dualCheck=0.0;
    16511671      if (totalThru+thruThis>=dualCheck) {
     
    24722492  int numberXColumns = info->numberXColumns();
    24732493  int numberXRows = info->numberXRows();
    2474   int iSequence;
    24752494  int start=numberXColumns+numberXRows;
    24762495  numberDualInfeasibilities_=0;
    24772496  sumDualInfeasibilities_=0.0;
    24782497  const int * which = info->quadraticSequence();
    2479   for (iSequence=0;iSequence<numberXColumns;iSequence++) {
    2480     int jSequence = which[iSequence];
    2481     double value=dj_[iSequence];
    2482     ClpSimplex::Status status = getColumnStatus(iSequence);
     2498  // Compute objective function from scratch
     2499  const CoinPackedMatrix * quadratic =
     2500    info->originalModel()->quadraticObjective();
     2501  const int * columnQuadratic = quadratic->getIndices();
     2502  const int * columnQuadraticStart = quadratic->getVectorStarts();
     2503  const int * columnQuadraticLength = quadratic->getVectorLengths();
     2504  const double * quadraticElement = quadratic->getElements();
     2505  const double * originalCost = info->originalObjective();
     2506  int iColumn;
     2507  objectiveValue_=0.0;
     2508  double infeasCost=0.0;
     2509  double linearCost=0.0;
     2510  for (iColumn=0;iColumn<numberXColumns;iColumn++) {
     2511    double valueI = solution_[iColumn];
     2512    linearCost += valueI*originalCost[iColumn];
     2513    double diff =cost_[iColumn]-originalCost[iColumn];
     2514    // WE are always feasible so look at low not up!
     2515    if (diff>0.0) {
     2516      double above = valueI-lower_[iColumn];
     2517      assert(above>=0.0);
     2518      infeasCost += diff*above;
     2519    } else if (diff<0.0) {
     2520      double below = upper_[iColumn]-valueI;
     2521      assert(below>=0.0);
     2522      infeasCost -= diff*below;
     2523    }
     2524    int j;
     2525    for (j=columnQuadraticStart[iColumn];
     2526         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     2527      int jColumn = columnQuadratic[j];
     2528      double valueJ = solution_[jColumn];
     2529      double elementValue = 0.5*quadraticElement[j];
     2530      objectiveValue_ += valueI*valueJ*elementValue;
     2531    }
     2532    int jSequence = which[iColumn];
     2533    double value=dj_[iColumn];
     2534    ClpSimplex::Status status = getColumnStatus(iColumn);
    24832535   
    24842536    switch(status) {
     
    24892541        if (getColumnStatus(jSequence)==basic)
    24902542          handler_->message(CLP_QUADRATIC_BOTH,messages_)
    2491             <<"Structural"<<iSequence
    2492             <<solution_[iSequence]<<jSequence<<solution_[jSequence]
     2543            <<"Structural"<<iColumn
     2544            <<solution_[iColumn]<<jSequence<<solution_[jSequence]
    24932545            <<CoinMessageEol;
    24942546      }
     
    25172569  }
    25182570  int offset = numberXColumns;
    2519   for (iSequence=0;iSequence<numberXRows;iSequence++) {
    2520     int jSequence = iSequence+offset;
    2521     double value=dj_[iSequence+numberColumns_];
    2522     ClpSimplex::Status status = getRowStatus(iSequence);
     2571  int iRow;
     2572  for (iRow=0;iRow<numberXRows;iRow++) {
     2573    int iColumn = iRow + numberColumns_;
     2574    double diff =cost_[iColumn];
     2575    double valueI = solution_[iColumn];
     2576    if (diff>0.0) {
     2577      double above = valueI-lower_[iColumn];
     2578      assert(above>=0.0);
     2579      infeasCost += diff*above;
     2580    } else if (diff<0.0) {
     2581      double below = upper_[iColumn]-valueI;
     2582      assert(below>=0.0);
     2583      infeasCost -= diff*below;
     2584    }
     2585    int jSequence = iRow+offset;
     2586    double value=dj_[iRow+numberColumns_];
     2587    ClpSimplex::Status status = getRowStatus(iRow);
    25232588   
    25242589    switch(status) {
     
    25272592      if (getColumnStatus(jSequence)==basic)
    25282593        printf("Row %d (%g) and %d (%g) both basic\n",
    2529                iSequence,solution_[iSequence],jSequence,solution_[jSequence]);
     2594               iRow,solution_[iColumn],jSequence,solution_[jSequence]);
    25302595      break;
    25312596    case ClpSimplex::isFixed:
     
    25512616    }
    25522617  }
     2618  objectiveValue_ += linearCost+infeasCost;
     2619  assert (infeasCost>=0.0);
    25532620  sumOfRelaxedDualInfeasibilities_ =sumDualInfeasibilities_;
     2621  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_;
    25542622  return numberDualInfeasibilities_;
    25552623}
     
    29473015      setColumnStatus(iColumn,status);
    29483016      if (status==basic) {
    2949         assert(fabs(value)<1.0e-7);
     3017        assert(fabs(value)<1.0e-2);
    29503018      }
    29513019      int j;
     
    29963064    validPhase_(-1),
    29973065    validSolution_(NULL),
     3066    djWeight_(NULL),
    29983067    numberXRows_(-1),
    29993068    numberXColumns_(-1),
     
    30143083    validPhase_(-1),
    30153084    validSolution_(NULL),
     3085    djWeight_(NULL),
    30163086    numberXRows_(-1),
    30173087    numberXColumns_(-1),
     
    30293099      backSequence_[i]=i;
    30303100    }
     3101    int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     3102    int numberRows = numberXRows_+numberQuadraticColumns_;
     3103    int size = numberRows+numberColumns;
     3104    djWeight_ = new double [size];
     3105    for (i=0;i<size;i++)
     3106      djWeight_[i]=1.0;
    30313107  }
    30323108}
     
    30383114  delete [] currentSolution_;
    30393115  delete [] validSolution_;
     3116  delete [] djWeight_;
    30403117}
    30413118// Copy
     
    30763153      validSolution_ = NULL;
    30773154    }
     3155    if (rhs.djWeight_) {
     3156      djWeight_ = new double [size];
     3157      memcpy(djWeight_,rhs.djWeight_,size*sizeof(double));
     3158    } else {
     3159      djWeight_ = NULL;
     3160    }
    30783161  }
    30793162}
     
    30883171    delete [] currentSolution_;
    30893172    delete [] validSolution_;
     3173    delete [] djWeight_;
    30903174    currentSequenceIn_ = rhs.currentSequenceIn_;
    30913175    crucialSj_ = rhs.crucialSj_;
     
    31223206    } else {
    31233207      validSolution_ = NULL;
     3208    }
     3209    if (rhs.djWeight_) {
     3210      djWeight_ = new double [size];
     3211      memcpy(djWeight_,rhs.djWeight_,size*sizeof(double));
     3212    } else {
     3213      djWeight_ = NULL;
    31243214    }
    31253215  }
  • branches/pre/include/ClpSimplexPrimalQuadratic.hpp

    r186 r191  
    174174  /// Restore previous
    175175  void restoreStatus();
     176  ///Dj weights
     177  inline double * djWeight() const
     178  { return djWeight_;};
    176179  //@}
    177180   
     
    201204  /// Valid saved solution
    202205  double * validSolution_;
     206  /// Dj weights to stop looping
     207  double * djWeight_;
    203208  /// Number of original rows
    204209  int numberXRows_;
Note: See TracChangeset for help on using the changeset viewer.