Changeset 186 for branches/pre


Ignore:
Timestamp:
Jul 28, 2003 9:43:58 AM (18 years ago)
Author:
forrest
Message:

stuff

Location:
branches/pre
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r185 r186  
    691691                      ClpQuadraticInfo * info)
    692692{
    693   checkComplementarity (info);
     693  checkComplementarity (info,rowArray_[0],rowArray_[1]);
    694694  int crucialSj = info->crucialSj();
    695695  int returnCode=-1;
    … …  
    767767      sequenceIn_ = nextSequenceIn;
    768768      cleanupIteration=false;
     769      lowerIn_=lower_[sequenceIn_];
     770      upperIn_=upper_[sequenceIn_];
    769771      if (sequenceIn_<numberXColumns) {
    770772        int jSequence = which[sequenceIn_];
    … …  
    887889        if (pivotRow_>=0) {
    888890          // if stable replace in basis
    889           int updateStatus = factorization_->replaceColumn(rowArray_[2],
    890                                                            pivotRow_,
    891                                                            alpha_);
     891          int updateStatus = 0;
     892          if (result<20)
     893            updateStatus=factorization_->replaceColumn(rowArray_[2],
     894                                                       pivotRow_,
     895                                                       alpha_);
    892896          if (result>=10) {
    893897            updateStatus = max(updateStatus,result/10);
    894898            result = result%10;
    895             if (updateStatus>1)
     899            if (updateStatus>1) {
    896900              alpha_=0.0; // force re-factorization
     901              info->setSequenceIn(sequenceIn_);
     902            }
    897903          }
    898904          // if no pivots, bad update but reasonable alpha - take and invert
    … …  
    10391045          setStatus(sequenceOut_,isFree);
    10401046        }
    1041         checkComplementarity (info);
     1047        checkComplementarity (info,rowArray_[2],rowArray_[3]);
    10421048
    10431049        if (whatNext==1) {
    … …  
    11971203  const int * columnQuadraticLength = quadratic->getVectorLengths();
    11981204  const double * quadraticElement = quadratic->getElements();
    1199   const double * originalCost = info->originalObjective();
     1205  //const double * originalCost = info->originalObjective();
    12001206  // Use rhsArray
    12011207  rhsArray->clear();
    … …  
    12071213  int iSjRow=-1;
    12081214  double tj = 0.0;
     1215  // Change in objective will be theta*coeff1 + theta*theta*coeff2
     1216  double coeff1 = 0.0;
     1217  double coeff2 = 0.0;
     1218  double coeffSlack=0.0;
    12091219  for (iIndex=0;iIndex<number;iIndex++) {
    12101220    int iRow = which[iIndex];
    … …  
    12141224      index2[number2++]=iPivot;
    12151225      rhs[iPivot]=alpha;
     1226      coeff1 += alpha*cost_[iPivot];
    12161227      //printf("col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
    12171228    } else {
     1229      coeffSlack += alpha*cost_[iPivot];
    12181230      //printf("new col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
    12191231      if (iPivot==crucialSj) {
    … …  
    12281240    rhs[sequenceIn_]=way;
    12291241    printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     1242    coeff1 += way*cost_[sequenceIn_];
    12301243  } else {
    12311244    printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     1245    coeffSlack += way*cost_[sequenceIn_];
    12321246  }
    12331247  rhsArray->setNumElements(number2);
    1234   // Change in objective will be theta*coeff1 + theta*theta*coeff2
    1235   double coeff1 = 0.0;
    1236   double coeff2 = 0.0;
    12371248  if (numberIterations_>=0||cleanupIteration) {
    12381249    for (iIndex=0;iIndex<number2;iIndex++) {
    … …  
    12401251      //double valueI = solution_[iColumn];
    12411252      double alphaI = rhs[iColumn];
    1242       coeff1 += alphaI*originalCost[iColumn];
    12431253      int j;
    12441254      for (j=columnQuadraticStart[iColumn];
    … …  
    12791289      //rhs[iColumn]=0.0;
    12801290      printf("Column %d, alpha %g sol %g cost %g, contr %g\n",
    1281              iColumn,alphaI,valueI,originalCost[iColumn],
    1282              alphaI*originalCost[iColumn]);
    1283       coeff1 += alphaI*originalCost[iColumn];
     1291             iColumn,alphaI,valueI,cost_[iColumn],
     1292             alphaI*cost_[iColumn]);
    12841293      int j;
    12851294      for (j=columnQuadraticStart[iColumn];
    … …  
    12981307  }
    12991308  coeff2 *= 0.5;
    1300   printf("coefficients %g %g - dualIn %g\n",coeff1,coeff2,dualIn_);
     1309  printf("coefficients %g %g %g - dualIn %g\n",coeff1,coeff2,coeffSlack,dualIn_);
     1310  if (coeffSlack)
     1311    printf("trouble\n");
     1312  coeff1 += coeffSlack;
     1313  //coeff1 = way*dualIn_;
    13011314  int accuracyFlag=0;
    13021315  if (!cleanupIteration) {
     1316    assert (fabs(way*coeff1-dualIn_)<1.0e-4*(1.0+fabs(dualIn_)));
     1317    assert (way*coeff1*dualIn_>0.0);
    13031318    if (way*coeff1*dualIn_<0.0) {
    13041319      // bad
    … …  
    13151330      solution_[crucialSj]=dualIn_;
    13161331    }
     1332  } else if (dualIn_<0.0&&way*coeff1>1.0e-2) {
     1333    printf("bad pair\n");
     1334    accuracyFlag=1;
     1335  } else if (dualIn_>0.0&&way*coeff1<-1.0e-2) {
     1336    accuracyFlag=1;
     1337    printf("bad pair\n");
    13171338  }
    13181339  // interesting places are where derivative zero or sj goes to zero
    … …  
    13411362    printf("linear variable\n");
    13421363  }
    1343   maximumMovement = min(maximumMovement,d1);
    13441364  maximumMovement = min(maximumMovement,d2);
    1345   d2 = min(d1,d2);
     1365  if (d1>=0.0) {
     1366    maximumMovement = min(maximumMovement,d1);
     1367    d2 = min(d1,d2);
     1368  }
    13461369   
    13471370  rhsArray->clear();
    … …  
    14381461      // but make a bit more pessimistic
    14391462      dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
     1463      //dualCheck=0.0;
    14401464      if (totalThru+thruThis>=dualCheck) {
    14411465        // We should be pivoting in this batch
    … …  
    15161540          } else {
    15171541            dualCheck = - 2.0*coeff2*theta_ - coeff1-9999999;
     1542            //dualCheck =0.0;
    15181543            if (totalThru>=dualCheck) {
    15191544              break; // no point trying another loop
    … …  
    16571682    } else {
    16581683      // need to do something
    1659       abort();
     1684      accuracyFlag=2;
    16601685    }
    16611686  }
    … …  
    16991724    printf("Objective value %g\n",objectiveValue());
    17001725  }
    1701   return result+10*accuracyFlag;
     1726  if (accuracyFlag<2) {
     1727    return result+10*accuracyFlag;
     1728  } else {
     1729    pivotRow_=1; // make sure correct path
     1730    return 20;
     1731  }
    17021732}
    17031733/* Checks if finished.  Updates status */
    … …  
    17601790  createRim(4);
    17611791  gutsOfSolution(NULL,NULL);
    1762   checkComplementarity(info);
     1792  checkComplementarity(info,rowArray_[2],rowArray_[3]);
    17631793  // Double check reduced costs if no action
    17641794  if (progress->lastIterationNumber(0)==numberIterations_) {
    … …  
    17801810    // something may have changed
    17811811    gutsOfSolution(NULL,NULL);
    1782     checkComplementarity(info);
     1812    checkComplementarity(info,rowArray_[2],rowArray_[3]);
    17831813  }
    17841814  // really for free variables in
    … …  
    18241854      nonLinearCost_->checkInfeasibilities(primalTolerance_);
    18251855      gutsOfSolution(NULL,NULL);
    1826       checkComplementarity(info);
     1856      checkComplementarity(info,rowArray_[2],rowArray_[3]);
    18271857
    18281858      infeasibilityCost_=1.0e100;
    … …  
    18421872          cost_[i] *= 1.0e-100;
    18431873        gutsOfSolution(NULL,NULL);
    1844         checkComplementarity(info);
    18451874        nonLinearCost_=nonLinear;
    18461875        infeasibilityCost_=saveWeight;
    … …  
    18611890          nonLinearCost_->checkInfeasibilities(primalTolerance_);
    18621891          gutsOfSolution(NULL,NULL);
    1863           checkComplementarity(info);
     1892          checkComplementarity(info,rowArray_[2],rowArray_[3]);
    18641893          // so will exit
    18651894          infeasibilityCost_=1.0e30;
    … …  
    18791908          nonLinearCost_->checkInfeasibilities(0.0);
    18801909          gutsOfSolution(NULL,NULL);
    1881           checkComplementarity(info);
     1910          checkComplementarity(info,rowArray_[2],rowArray_[3]);
    18821911          problemStatus_=-1; //continue
    18831912        } else {
    … …  
    19131942          nonLinearCost_->checkInfeasibilities(oldTolerance);
    19141943          gutsOfSolution(NULL,NULL);
    1915           checkComplementarity(info);
     1944          checkComplementarity(info,rowArray_[2],rowArray_[3]);
    19161945          problemStatus_ = -1;
    19171946        } else {
    … …  
    19451974          createRim(4);
    19461975          gutsOfSolution(NULL,NULL);
    1947           checkComplementarity(info);
     1976          checkComplementarity(info,rowArray_[2],rowArray_[3]);
    19481977          problemStatus_=-1; //continue
    19491978        } else {
    … …  
    19942023// Fills in reduced costs
    19952024void
    1996 ClpSimplexPrimalQuadratic::createDjs (const ClpQuadraticInfo * info)
     2025ClpSimplexPrimalQuadratic::createDjs (const ClpQuadraticInfo * info,
     2026                            CoinIndexedVector * array1, CoinIndexedVector * array2)
    19972027{
    19982028  const int * which = info->quadraticSequence();
    … …  
    20072037  int iSequence;
    20082038  int start=numberXColumns+numberXRows;
     2039  int jSequence;
     2040  // Actually only need to do this after invert (use extra parameter to createDjs)
     2041  // or when infeasibilities change
     2042  // recode to do this later and update rather than recompute
     2043  // When it is "change" then we don't need b (original rhs)
     2044  if (nonLinearCost_->numberInfeasibilities()<-1) {
     2045  } else {
     2046    array1->checkClear();
     2047    array2->checkClear();
     2048    // update costs
     2049    double * storedCost = lower_+numberColumns_+numberXRows;
     2050    double * storedUpper = upper_ + numberColumns_+numberXRows;
     2051    double * storedValue = solution_ + numberColumns_+numberXRows;
     2052    int * index = array1->getIndices();
     2053    double * modifiedCost = array1->denseVector();
     2054    // zero out basic values
     2055    int iRow;
     2056    double * save = new double [numberRows_];
     2057    for (iRow=0;iRow<numberRows_;iRow++) {
     2058      int iPivot=pivotVariable_[iRow];
     2059      save[iRow]=solution_[iPivot];
     2060      solution_[iPivot] = 0.0;
     2061    }
     2062    // And make sure pi is zero if slack basic
     2063    for (iRow=0;iRow<numberXRows;iRow++) {
     2064      //if (getRowStatus(iRow)==basic)
     2065      //assert(!solution_[iRow+numberXColumns]);
     2066      if (getRowStatus(iRow)==basic)
     2067        solution_[iRow+numberXColumns]=0.0;
     2068    }
     2069    times(-1.0,columnActivityWork_,modifiedCost);
     2070   
     2071    int number=0;
     2072    // Now costs
     2073    for (iSequence=0;iSequence<numberXColumns;iSequence++) {
     2074      iRow = iSequence+numberXRows;
     2075      double value=cost_[iSequence];
     2076      if (fabs(value)>1.0e5) {
     2077        assert (getColumnStatus(iSequence)==basic);
     2078      }
     2079      storedCost[iSequence]=value;
     2080      storedUpper[iSequence]=value;
     2081      storedValue[iSequence]=value;
     2082      value += modifiedCost[iRow];
     2083      if (value) {
     2084        modifiedCost[iRow]=value;
     2085        index[number++]=iRow;
     2086      } else {
     2087        modifiedCost[iRow]=0.0;
     2088      }
     2089    }
     2090    for (iRow=0;iRow<numberXRows;iRow++) {
     2091      double value = modifiedCost[iRow]+rowActivityWork_[iRow];;
     2092      if (value) {
     2093        modifiedCost[iRow]=value;
     2094        index[number++]=iRow;
     2095      } else {
     2096        modifiedCost[iRow]=0.0;
     2097      }
     2098    }
     2099    array1->setNumElements(number);
     2100    // And slacks
     2101    // sign? - or does any of this make sense
     2102    for (iRow=0;iRow<numberXRows;iRow++) {
     2103      int iSequence  = iRow + numberColumns_;
     2104      double value = cost_[iSequence];
     2105      if (value) {
     2106        assert(getRowStatus(iRow)==basic);
     2107        // add in pi column
     2108        matrix_->add(this,array1,iRow+numberXColumns,value);
     2109      }
     2110    }
     2111    // create pricing vector
     2112    int k;
     2113    factorization_->updateColumn(array2,array1);
     2114    int n=array1->getNumElements();
     2115    for (k=0;k<n;k++) {
     2116      int iRow = index[k];
     2117      double value=modifiedCost[iRow];
     2118      int iSequence = pivotVariable_[iRow];
     2119      if (iSequence>=numberXColumns&&iSequence<numberColumns_)
     2120        solution_[iSequence] = value;
     2121      else
     2122        solution_[iSequence] = save[iRow];
     2123      //solution_[iSequence] = value;
     2124     
     2125    }
     2126    delete [] save;
     2127    array1->clear();
     2128    // Now modify pi values by slack costs to make djs
     2129    // Later save matrix_->add results
     2130    for (iRow=0;iRow<numberXRows;iRow++) {
     2131      int iSequence  = iRow + numberColumns_;
     2132      int iPi = iRow+numberXColumns;
     2133      solution_[iPi]-=cost_[iSequence];
     2134    }
     2135    // check looks okay
     2136    matrix_->times(1.0,solution_,modifiedCost);
     2137    // Back to pi
     2138    for (iRow=0;iRow<numberXRows;iRow++) {
     2139      int iSequence  = iRow + numberColumns_;
     2140      int iPi = iRow+numberXColumns;
     2141      solution_[iPi]+=cost_[iSequence];
     2142    }
     2143    double * rowSolution = solution_+numberColumns_;
     2144    for (k=0;k<numberRows_;k++) {
     2145      if (k<numberXRows) {
     2146        if (fabs(modifiedCost[k]-rowSolution[k])>1.0e-3*(1.0+fabs(rowSolution[k])))
     2147          printf("bad row %d %g %g\n",k,modifiedCost[k],rowSolution[k]);
     2148      } else {
     2149        rowSolution[k]=modifiedCost[k];
     2150        lower_[k+numberColumns_]=modifiedCost[k];
     2151        upper_[k+numberColumns_]=modifiedCost[k];
     2152      }
     2153      modifiedCost[k]=0.0;
     2154    }
     2155  }
    20092156  for (iSequence=0;iSequence<numberXColumns;iSequence++) {
    20102157    int jSequence = which[iSequence];
    … …  
    20142161      value = solution_[jSequence];
    20152162    } else {
    2016       value=cost_[iSequence];
     2163      value=0.0;
     2164      //value=cost_[iSequence];
    20172165      int j;
    20182166      for (j=columnStart[iSequence];j<columnStart[iSequence]+columnLength[iSequence]; j++) {
    … …  
    20252173  // And slacks
    20262174  // Value of sj is - value of pi
    2027   int firstSlack = numberColumns_;
    2028   int jSequence;
    20292175  for (jSequence=0;jSequence<numberXRows;jSequence++) {
    2030     int iSequence  = jSequence + firstSlack;
     2176    int iSequence  = jSequence + numberColumns_;
    20312177    int iPi = jSequence+numberXColumns;
    20322178    double value = solution_[iPi];
    2033     dj_[iSequence]=value;
     2179    dj_[iSequence]=value+cost_[iSequence];
    20342180  }
    20352181}
    20362182
    20372183int
    2038 ClpSimplexPrimalQuadratic::checkComplementarity (const
    2039                                                ClpQuadraticInfo * info)
     2184ClpSimplexPrimalQuadratic::checkComplementarity (const  ClpQuadraticInfo * info,
     2185                            CoinIndexedVector * array1, CoinIndexedVector * array2)
    20402186{
    2041   createDjs(info);
     2187  createDjs(info,array1,array2);
    20422188  int numberXColumns = info->numberXColumns();
    20432189  int numberXRows = info->numberXRows();
    … …  
    20932239     
    20942240    case ClpSimplex::basic:
    2095       if (getRowStatus(jSequence)==basic)
     2241      if (getColumnStatus(jSequence)==basic)
    20962242        printf("Row %d (%g) and %d (%g) both basic\n",
    20972243               iSequence,solution_[iSequence],jSequence,solution_[jSequence]);
    … …  
    22332379    columnUpper2[iColumn]=columnUpper[iColumn];
    22342380    solution2[iColumn]=solution[iColumn];
     2381    // Put in cost so we know it
     2382    objective2[iColumn]=objective[iColumn];
    22352383    int j;
    22362384    for (j=columnStart[iColumn];
    … …  
    24512599  }
    24522600  printf("Objective value %g\n",objValue);
    2453   for (iColumn=0;iColumn<newNumberColumns;iColumn++)
    2454     printf("%d %g\n",iColumn,solution2[iColumn]);
     2601  //for (iColumn=0;iColumn<newNumberColumns;iColumn++)
     2602  //printf("%d %g\n",iColumn,solution2[iColumn]);
    24552603  return (ClpSimplexPrimalQuadratic *) model2;
    24562604}

  • branches/pre/Test/unitTest.cpp

    r183 r186  
    1717#include "ClpFactorization.hpp"
    1818#include "ClpSimplex.hpp"
     19#include "ClpLinearObjective.hpp"
    1920#include "ClpDualRowSteepest.hpp"
    2021#include "ClpDualRowDantzig.hpp"
    … …  
    11181119      start[i+1]=i+1;
    11191120      column[i]=i;
    1120       element[i]=1.0e-6;
     1121      element[i]=1.0e-1;
    11211122      element[i]=0.0;
    11221123    }
    … …  
    11421143    //fn = "share2qpb";
    11431144    m.readMps(fn.c_str(),"mps");
    1144     ClpSimplex solution;
    1145     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
     1145    ClpSimplex model;
     1146    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    11461147                         m.getObjCoefficients(),
    11471148                         m.getRowLower(),m.getRowUpper());
    1148     solution.dual();
     1149    model.dual();
    11491150    // get quadratic part
    11501151    int * start=NULL;
    … …  
    11531154    m.readQuadraticMps(NULL,start,column,element,2);
    11541155    // Load up objective
    1155     solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
     1156    model.loadQuadraticObjective(model.numberColumns(),start,column,element);
    11561157    delete [] start;
    11571158    delete [] column;
    11581159    delete [] element;
    1159     solution.quadraticSLP(50,1.0e-4);
    1160     solution.setLogLevel(63);
    1161     solution.quadraticPrimal(1);
    1162     double objValue = solution.getObjValue();
     1160#if 0
     1161    model.quadraticSLP(50,1.0e-4);
     1162#else
     1163    // Get feasible
     1164    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
     1165    int numberColumns=model.numberColumns();
     1166    ClpLinearObjective zeroObjective(NULL,numberColumns);
     1167    model.setObjective(&zeroObjective);
     1168    model.dual();
     1169    model.setObjective(saveObjective);
     1170    delete saveObjective;
     1171#endif
     1172    model.setLogLevel(63);
     1173    model.quadraticPrimal(1);
     1174    double objValue = model.getObjValue();
    11631175    CoinRelFltEq eq(1.0e-4);
    11641176    assert(eq(objValue,-400.92));

  • branches/pre/include/ClpFactorization.hpp

    r180 r186  
    8282                     CoinIndexedVector * regionSparse2,
    8383                     bool noPermute=false) const;
    84   /** Updates one column transpose (BTRAN)
    85       ** For large problems you should ALWAYS know where the nonzeros
    86       are, so please try and migrate to previous method after you
    87       have got code working using this simple method - thank you!
    88       (the only exception is if you know input is dense e.g. dense objective)
    89       returns number of nonzeros */
    90   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    91                                  double array[] ) const;
    9284  /** Updates one column (BTRAN) from region2
    9385      region1 starts as zero and is zero at end */

  • branches/pre/include/ClpSimplexPrimalQuadratic.hpp

    r185 r186  
    6060                   ClpQuadraticInfo & info);
    6161  /// Checks complementarity and computes infeasibilities
    62   int checkComplementarity (const ClpQuadraticInfo * info);
     62  int checkComplementarity (const ClpQuadraticInfo * info,
     63                            CoinIndexedVector * array1, CoinIndexedVector * array2);
    6364  /// Fills in reduced costs
    64   void createDjs (const ClpQuadraticInfo * info);
     65  void createDjs (const ClpQuadraticInfo * info,
     66                  CoinIndexedVector * array1, CoinIndexedVector * array2);
    6567  /** Main part.
    6668      phase - 0 normal, 1 getting complementary solution,
Note: See TracChangeset for help on using the changeset viewer.