Changeset 678 for branches/devel


Ignore:
Timestamp:
Jul 9, 2007 4:52:34 PM (12 years ago)
Author:
forrest
Message:

for tighter linear constraints in bilinear

Location:
branches/devel/Cbc/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Cbc/src/CbcLinked.cpp

    r667 r678  
    181181  //sprintf(temp,"cc%d",iPass);
    182182  //writeMps(temp);
     183  //writeMps("tight");
     184  //exit(33);
    183185  //printf("wrote cc%d\n",iPass);
    184186  OsiClpSolverInterface::initialSolve();
     
    14971499    }
    14981500  }
     1501  delete [] which;
     1502#if 1
     1503  addTighterConstraints();
     1504#endif
     1505}
     1506// Add reformulated bilinear constraints
     1507void
     1508OsiSolverLink::addTighterConstraints()
     1509{
     1510  // This is first attempt - for now get working on trimloss
     1511  int numberW=0;
     1512  int * xW = new int[numberObjects_];
     1513  int * yW = new int[numberObjects_];
     1514  // Points to firstlambda
     1515  int * wW = new int[numberObjects_];
     1516  // Coefficient
     1517  double * alphaW = new double[numberObjects_];
     1518  // Objects
     1519  OsiBiLinear ** objW = new OsiBiLinear * [numberObjects_];
     1520  int numberColumns = getNumCols();
     1521  int firstLambda=numberColumns;
     1522  // set up list (better to rethink and do properly as column ordered)
     1523  int * list = new int[numberColumns];
     1524  memset(list,0,numberColumns*sizeof(int));
     1525  int i;
     1526  for ( i =0;i<numberObjects_;i++) {
     1527    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     1528    if (obj) {
     1529      //obj->setBranchingStrategy(4); // ***** temp
     1530      objW[numberW]=obj;
     1531      xW[numberW]=obj->xColumn();
     1532      yW[numberW]=obj->yColumn();
     1533      list[xW[numberW]]=1;
     1534      list[yW[numberW]]=1;
     1535      wW[numberW]=obj->firstLambda();
     1536      firstLambda = CoinMin(firstLambda,obj->firstLambda());
     1537      alphaW[numberW]=obj->coefficient();
     1538      //assert (alphaW[numberW]==1.0); // fix when occurs
     1539      numberW++;
     1540    }
     1541  }
     1542  int nList = 0;
     1543  for (i=0;i<numberColumns;i++) {
     1544    if (list[i])
     1545      list[nList++]=i;
     1546  }
     1547  // set up mark array
     1548  char * mark = new char [firstLambda*firstLambda];
     1549  memset(mark,0,firstLambda*firstLambda);
     1550  for (i=0;i<numberW;i++) {
     1551    int x = xW[i];
     1552    int y = yW[i];
     1553    mark[x*firstLambda+y]=1;
     1554    mark[y*firstLambda+x]=1;
     1555  }
     1556  int numberRows2 = originalRowCopy_->getNumRows();
     1557  int * addColumn = new int [numberColumns];
     1558  double * addElement = new double [numberColumns];
     1559  assert (objectiveRow_<0); // fix when occurs
     1560  for (int iRow=0;iRow<numberRows2;iRow++) {
     1561    for (int iList=0;iList<nList;iList++) {
     1562      int kColumn = list[iList];
     1563      const double * columnLower = getColLower();
     1564      //const double * columnUpper = getColUpper();
     1565      const double * rowLower = getRowLower();
     1566      const double * rowUpper = getRowUpper();
     1567      const CoinPackedMatrix * rowCopy = getMatrixByRow();
     1568      const double * element = rowCopy->getElements();
     1569      const int * column = rowCopy->getIndices();
     1570      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     1571      const int * rowLength = rowCopy->getVectorLengths();
     1572      CoinBigIndex j;
     1573      int numberElements = rowLength[iRow];
     1574      int n=0;
     1575      for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
     1576        int iColumn = column[j];
     1577        if (iColumn>=firstLambda) {
     1578          // no good
     1579          n=-1;
     1580          break;
     1581        }
     1582        if (mark[iColumn*firstLambda+kColumn])
     1583          n++;
     1584      }
     1585      if (n==numberElements) {
     1586        printf("can add row %d\n",iRow);
     1587        assert (columnLower[kColumn]>=0); // might be able to fix
     1588        int xColumn=kColumn;
     1589        n=0;
     1590        for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
     1591          int yColumn = column[j];
     1592          int k;
     1593          for (k=0;k<numberW;k++) {
     1594            if ((xW[k]==yColumn&&yW[k]==xColumn)||
     1595                (yW[k]==yColumn&&xW[k]==xColumn))
     1596              break;
     1597          }
     1598          assert (k<numberW);
     1599          if (xW[k]!=xColumn) {
     1600            int temp=xColumn;
     1601            xColumn=yColumn;
     1602            yColumn=temp;
     1603          }
     1604          int start = wW[k];
     1605          double value = element[j];
     1606          for (int kk=0;kk<4;kk++) {
     1607            // Dummy value
     1608            addElement[n]= 1.0e-19;
     1609            addColumn[n++]=start+kk;
     1610          }
     1611          // Tell object about this
     1612          objW[k]->addExtraRow(matrix_->getNumRows(),value);
     1613        }
     1614        addColumn[n++] = kColumn;
     1615        double lo = rowLower[iRow];
     1616        double up = rowUpper[iRow];
     1617        if (lo>-1.0e20) {
     1618          addElement[n-1]=-lo;
     1619          if (lo==up)
     1620            addRow(n,addColumn,addElement,0.0,0.0);
     1621          else
     1622            addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
     1623          matrix_->appendRow(n,addColumn,addElement);
     1624        }
     1625        if (up<1.0e20&&up>lo) {
     1626          addElement[n-1]=-up;
     1627          addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
     1628          matrix_->appendRow(n,addColumn,addElement);
     1629        }
     1630      }
     1631    }
     1632  }
    14991633#if 0
    1500   //int fake[10]={3,4,2,5,5,3,1,11,2,0};
    1501   int fake[10]={11,5,5,4,3,3,2,2,1,0};
    1502   for (int kk=0;kk<10;kk++) {
    1503     setColUpper(kk,fake[kk]);
    1504     setColLower(kk,fake[kk]);
     1634  // possibly do bounds
     1635  for (int iColumn=0;iColumn<firstLambda;iColumn++) {
     1636    for (int iList=0;iList<nList;iList++) {
     1637      int kColumn = list[iList];
     1638      const double * columnLower = getColLower();
     1639      const double * columnUpper = getColUpper();
     1640      if (mark[iColumn*firstLambda+kColumn]) {
     1641        printf("can add column %d\n",iColumn);
     1642        assert (columnLower[kColumn]>=0); // might be able to fix
     1643        int xColumn=kColumn;
     1644        int yColumn = iColumn;
     1645        int k;
     1646        for (k=0;k<numberW;k++) {
     1647          if ((xW[k]==yColumn&&yW[k]==xColumn)||
     1648              (yW[k]==yColumn&&xW[k]==xColumn))
     1649            break;
     1650        }
     1651        assert (k<numberW);
     1652        if (xW[k]!=xColumn) {
     1653          int temp=xColumn;
     1654          xColumn=yColumn;
     1655          yColumn=temp;
     1656        }
     1657        int start = wW[k];
     1658        int n=0;
     1659        for (int kk=0;kk<4;kk++) {
     1660          // Dummy value
     1661          addElement[n]= 1.0e-19;
     1662          addColumn[n++]=start+kk;
     1663        }
     1664        // Tell object about this
     1665        objW[k]->addExtraRow(matrix_->getNumRows(),1.0);
     1666        addColumn[n++] = kColumn;
     1667        double lo = columnLower[iColumn];
     1668        double up = columnUpper[iColumn];
     1669        if (lo>-1.0e20) {
     1670          addElement[n-1]=-lo;
     1671          if (lo==up)
     1672            addRow(n,addColumn,addElement,0.0,0.0);
     1673          else
     1674            addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
     1675          matrix_->appendRow(n,addColumn,addElement);
     1676        }
     1677        if (up<1.0e20&&up>lo) {
     1678          addElement[n-1]=-up;
     1679          addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
     1680          matrix_->appendRow(n,addColumn,addElement);
     1681        }
     1682      }
     1683    }
    15051684  }
    15061685#endif
    1507   delete [] which;
     1686  delete [] xW;
     1687  delete [] yW;
     1688  delete [] wW;
     1689  delete [] alphaW;
     1690  delete [] addColumn;
     1691  delete [] addElement;
     1692  delete [] mark;
     1693  delete [] list;
     1694  delete [] objW;
    15081695}
    15091696// Set all biLinear priorities on x-x variables
     
    44414628    xyRow_(-1),
    44424629    convexity_(-1),
     4630    numberExtraRows_(0),
     4631    multiplier_(NULL),
     4632    extraRow_(NULL),
    44434633    chosen_(-1)
    44444634{
     
    44694659    xyRow_(xyRow),
    44704660    convexity_(-1),
     4661    numberExtraRows_(0),
     4662    multiplier_(NULL),
     4663    extraRow_(NULL),
    44714664    chosen_(-1)
    44724665{
     
    46804873    xyRow_(xyRow),
    46814874    convexity_(-1),
     4875    numberExtraRows_(0),
     4876    multiplier_(NULL),
     4877    extraRow_(NULL),
    46824878    chosen_(-1)
    46834879{
     
    48965092   xyRow_(rhs.xyRow_),
    48975093   convexity_(rhs.convexity_),
     5094   numberExtraRows_(rhs.numberExtraRows_),
     5095   multiplier_(NULL),
     5096   extraRow_(NULL),
    48985097   chosen_(rhs.chosen_)
    48995098{
     5099  if (numberExtraRows_) {
     5100    multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
     5101    extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
     5102  }
    49005103}
    49015104
     
    49315134    xyRow_ = rhs.xyRow_;
    49325135    convexity_ = rhs.convexity_;
     5136    numberExtraRows_ = rhs.numberExtraRows_;
     5137    delete [] multiplier_;
     5138    delete [] extraRow_;
     5139    if (numberExtraRows_) {
     5140      multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
     5141      extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
     5142    } else {
     5143      multiplier_ = NULL;
     5144      extraRow_ = NULL;
     5145    }
    49335146    chosen_ = rhs.chosen_;
    49345147  }
     
    49395152OsiBiLinear::~OsiBiLinear ()
    49405153{
     5154  delete [] multiplier_;
     5155  delete [] extraRow_;
     5156}
     5157// Adds in data for extra row with variable coefficients
     5158void
     5159OsiBiLinear::addExtraRow(int row, double multiplier)
     5160{
     5161  int * tempI = new int [numberExtraRows_+1];
     5162  double * tempD = new double [numberExtraRows_+1];
     5163  memcpy(tempI,extraRow_,numberExtraRows_*sizeof(int));
     5164  memcpy(tempD,multiplier_,numberExtraRows_*sizeof(double));
     5165  tempI[numberExtraRows_]=row;
     5166  tempD[numberExtraRows_]=multiplier;
     5167  numberExtraRows_++;
     5168  delete [] extraRow_;
     5169  extraRow_ = tempI;
     5170  delete [] multiplier_;
     5171  multiplier_ = tempD;
    49415172}
    49425173static bool testCoarse=true;
     
    51215352    }
    51225353    xyLambda /= coefficient_;
     5354  }
     5355  if (0) {
     5356    // only true with positive values
     5357    // see if all convexification constraints OK with true
     5358    assert (xyTrue+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]);
     5359    assert (xyTrue+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]);
     5360    assert (xyTrue-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]);
     5361    assert (xyTrue-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]);
     5362    // see if all convexification constraints OK with lambda version
     5363#if 1
     5364    assert (xyLambda+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]);
     5365    assert (xyLambda+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]);
     5366    assert (xyLambda-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]);
     5367    assert (xyLambda-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]);
     5368#endif
     5369    // see if other bound stuff true
     5370    assert (xyLambda+1.0e-5>xB[0]*y);
     5371    assert (xyLambda+1.0e-5>yB[0]*x);
     5372    assert (xyLambda-1.0e-5<xB[1]*y);
     5373    assert (xyLambda-1.0e-5<yB[1]*x);
     5374#define SIZE 2
     5375    if (yColumn_==xColumn_+SIZE) {
     5376#if SIZE==6
     5377      double bMax = 2200.0;
     5378      double bMin = bMax - 100.0;
     5379      double b[] = {330.0,360.0,380.0,430.0,490.0,530.0};
     5380#elif SIZE==2
     5381      double bMax = 1900.0;
     5382      double bMin = bMax - 200.0;
     5383      double b[] = {460.0,570.0};
     5384#else
     5385      abort();
     5386#endif
     5387      double sum =0.0;
     5388      double sum2 =0.0;
     5389      int m=xColumn_;
     5390      double x = info->solution_[m];
     5391      double xB[2];
     5392      double yB[2];
     5393      xB[0]=info->lower_[m];
     5394      xB[1]=info->upper_[m];
     5395      for (int i=0;i<SIZE*SIZE;i+=SIZE) {
     5396        int n = i+SIZE+m;
     5397        double y = info->solution_[n];
     5398        yB[0]=info->lower_[n];
     5399        yB[1]=info->upper_[n];
     5400        int firstLambda=SIZE*SIZE+2*SIZE+4*i+4*m;
     5401        double xyLambda=0.0;
     5402        for (int j=0;j<4;j++) {
     5403          int iX = j>>1;
     5404          int iY = j&1;
     5405          xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j];
     5406        }
     5407        sum += xyLambda*b[i/SIZE];
     5408        double xyTrue = x*y;
     5409        sum2 += xyTrue*b[i/SIZE];
     5410      }
     5411      if (sum>bMax*x+1.0e-5||sum<bMin*x-1.0e-5) {
     5412        //if (sum<bMax*x+1.0e-5&&sum>bMin*x-1.0e-5) {
     5413        printf("bmin*x %g b*w %g bmax*x %g (true) %g\n",bMin*x,sum,bMax*x,sum2);
     5414        printf("m %d lb %g value %g up %g\n",
     5415               m,xB[0],x,xB[1]);
     5416        sum=0.0;
     5417        for (int i=0;i<SIZE*SIZE;i+=SIZE) {
     5418          int n = i+SIZE+m;
     5419          double y = info->solution_[n];
     5420          yB[0]=info->lower_[n];
     5421          yB[1]=info->upper_[n];
     5422          printf("n %d lb %g value %g up %g\n",
     5423                 n,yB[0],y,yB[1]);
     5424          int firstLambda=SIZE*SIZE+2*SIZE+4*i+m*4;
     5425          double xyLambda=0.0;
     5426          for (int j=0;j<4;j++) {
     5427            int iX = j>>1;
     5428            int iY = j&1;
     5429            xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j];
     5430            printf("j %d l %d new xylambda %g ",j,firstLambda+j,xyLambda);
     5431          }
     5432          sum += xyLambda*b[i/SIZE];
     5433          printf(" - sum now %g\n",sum);
     5434        }
     5435      }
     5436      if (sum2>bMax*x+1.0e-5||sum2<bMin*x-1.0e-5) {
     5437        printf("bmin*x %g b*x*y %g bmax*x %g (estimate) %g\n",bMin*x,sum2,bMax*x,sum);
     5438        printf("m %d lb %g value %g up %g\n",
     5439               m,xB[0],x,xB[1]);
     5440        sum2=0.0;
     5441        for (int i=0;i<SIZE*SIZE;i+=SIZE) {
     5442          int n = i+SIZE+m;
     5443          double y = info->solution_[n];
     5444          yB[0]=info->lower_[n];
     5445          yB[1]=info->upper_[n];
     5446          printf("n %d lb %g value %g up %g\n",
     5447                 n,yB[0],y,yB[1]);
     5448          double xyTrue = x*y;
     5449          sum2 += xyTrue*b[i/SIZE];
     5450          printf("xyTrue %g - sum now %g\n",xyTrue,sum2);
     5451        }
     5452      }
     5453    }
    51235454  }
    51245455  // If pseudo shadow prices then see what would happen
     
    58416172  const int * row = matrix->getIndices();
    58426173  const CoinBigIndex * columnStart = matrix->getVectorStarts();
    5843   //const int * columnLength = matrix->getVectorLengths();
     6174  const int * columnLength = matrix->getVectorLengths();
    58446175  // order is LxLy, LxUy, UxLy and UxUy
    58456176  double xB[2];
     
    58626193    double y = yB[iY];
    58636194    CoinBigIndex k = columnStart[j+firstLambda_];
     6195    CoinBigIndex last = k+columnLength[j+firstLambda_];
    58646196    double value;
    58656197    // xy
     
    58966228      element[k++]=value;
    58976229      numberUpdated++;
     6230    }
     6231    // Do extra rows
     6232    for (int i=0;i<numberExtraRows_;i++) {
     6233      int iRow = extraRow_[i];
     6234      for (;k<last;k++) {
     6235        if (row[k]==iRow)
     6236          break;
     6237      }
     6238      assert (k<last);
     6239      element[k++] = x*y*multiplier_[i];
    58986240    }
    58996241  }
  • branches/devel/Cbc/src/CbcLinked.hpp

    r642 r678  
    115115  /// Analyze constraints to see which are convex (quadratic)
    116116  void analyzeObjects();
     117  /// Add reformulated bilinear constraints
     118  void addTighterConstraints();
    117119  /// Objective value of best solution found internally
    118120  inline double bestObjectiveValue() const
     
    842844  /// Compute lambdas (third entry in each .B is current value) (nonzero if bad)
    843845  double computeLambdas(const double xB[3], const double yB[3],const double xybar[4],double lambda[4]) const;
     846  /// Adds in data for extra row with variable coefficients
     847  void addExtraRow(int row, double multiplier);
    844848
    845849protected:
     
    899903  /// Convexity row
    900904  int convexity_;
     905  /// Number of extra rows (coefficients to be modified)
     906  int numberExtraRows_;
     907  /// Multiplier for coefficient on row
     908  double * multiplier_;
     909  /// Row number
     910  int * extraRow_;
    901911  /// Which chosen -1 none, 0 x, 1 y
    902912  mutable short chosen_;
  • branches/devel/Cbc/src/CbcSolver.cpp

    r674 r678  
    964964  } else {
    965965    delete finalModel;
     966    printf("can't make knapsacks - did you set fixedPriority (extra1)\n");
    966967    return NULL;
    967968  }
     
    32333234                        babModel->assignSolver(solver);
    32343235                        testOsiOptions=0;
     3236                        // allow gomory
     3237                        complicatedInteger=0;
    32353238                        // Priorities already done
    32363239                        free(info.priorities);
     
    39703973                        }
    39713974                      }
    3972                     }
    3973                     if (priorities) {
    3974                       int iPriority = priorities[iColumn];
    3975                       if (iPriority>0)
    3976                         objects[iObj]->setPriority(iPriority);
    3977                     }
    3978                     if (logLevel>2)
    3979                       printf("Obj %d is int? - priority %d\n",iObj,objects[iObj]->priority());
    3980                     if (pseudoUp&&pseudoUp[iColumn]) {
    3981                       abort();
     3975                      if (priorities) {
     3976                        int iPriority = priorities[iColumn];
     3977                        if (iPriority>0)
     3978                          objects[iObj]->setPriority(iPriority);
     3979                      }
     3980                      if (logLevel>2)
     3981                        printf("Obj %d is int? - priority %d\n",iObj,objects[iObj]->priority());
     3982                      if (pseudoUp&&pseudoUp[iColumn]) {
     3983                        abort();
     3984                      }
    39823985                    }
    39833986                  }
Note: See TracChangeset for help on using the changeset viewer.