Changeset 1120 for stable


Ignore:
Timestamp:
Sep 25, 2007 10:13:31 AM (12 years ago)
Author:
forrest
Message:

add find network code

Location:
stable/1.6/Clp
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • stable/1.6/Clp/src/ClpModel.cpp

    r1060 r1120  
    35813581  return coinModel;
    35823582}
     3583/* Find a network subset.
     3584   rotate array should be numberRows.  On output
     3585   -1 not in network
     3586   0 in network as is
     3587   1 in network with signs swapped
     3588  Returns number of network rows (positive if exact network, negative if needs extra row)
     3589  From Gulpinar, Gutin, Maros and Mitra
     3590*/
     3591int
     3592ClpModel::findNetwork(char * rotate,double fractionNeeded)
     3593{
     3594  int * mapping = new int [numberRows_];
     3595  // Get column copy
     3596  CoinPackedMatrix * columnCopy = matrix();
     3597  // Get a row copy in standard format
     3598  CoinPackedMatrix * copy = new CoinPackedMatrix();
     3599  copy->reverseOrderedCopyOf(*columnCopy);
     3600  // make sure ordered and no gaps
     3601  copy->cleanMatrix();
     3602  // get matrix data pointers
     3603  const int * columnIn = copy->getIndices();
     3604  const CoinBigIndex * rowStartIn = copy->getVectorStarts();
     3605  const int * rowLength = copy->getVectorLengths();
     3606  const double * elementByRowIn = copy->getElements();
     3607  int iRow,iColumn;
     3608  int numberEligible=0;
     3609  int numberIn=0;
     3610  int numberElements=0;
     3611  for (iRow=0;iRow<numberRows_;iRow++) {
     3612    bool possible=true;
     3613    mapping[iRow]=-1;
     3614    rotate[iRow]=-1;
     3615    for (CoinBigIndex j=rowStartIn[iRow];j<rowStartIn[iRow]+rowLength[iRow];j++) {
     3616      //int iColumn = column[j];
     3617      double value = elementByRowIn[j];
     3618      if (fabs(value)!=1.0) {
     3619        possible=false;
     3620        break;
     3621      }
     3622    }
     3623    if (rowLength[iRow]&&possible) {
     3624      mapping[iRow]=numberEligible;
     3625      numberEligible++;
     3626      numberElements+=rowLength[iRow];
     3627    }
     3628  }
     3629  if (numberEligible<fractionNeeded*numberRows_) {
     3630    delete [] mapping;
     3631    return 0;
     3632  }
     3633  // create arrays
     3634  int * eligible = new int [numberRows_];
     3635  int * column = new int [numberElements];
     3636  CoinBigIndex * rowStart = new CoinBigIndex [numberEligible+1];
     3637  char * elementByRow = new char [numberElements];
     3638  numberEligible=0;
     3639  numberElements=0;
     3640  rowStart[0]=0;
     3641  for (iRow=0;iRow<numberRows_;iRow++) {
     3642    if (mapping[iRow]<0)
     3643      continue;
     3644    assert (numberEligible==mapping[iRow]);
     3645    rotate[numberEligible]=0;
     3646    for (CoinBigIndex j=rowStartIn[iRow];j<rowStartIn[iRow]+rowLength[iRow];j++) {
     3647      column[numberElements] = columnIn[j];
     3648      double value = elementByRowIn[j];
     3649      if (value==1.0)
     3650        elementByRow[numberElements++]=1;
     3651      else
     3652        elementByRow[numberElements++]=-1;
     3653    }
     3654    numberEligible++;
     3655    rowStart[numberEligible]=numberElements;
     3656  }
     3657  // get rid of copy to save space
     3658  delete copy;
     3659  const int * rowIn = columnCopy->getIndices();
     3660  const CoinBigIndex * columnStartIn = columnCopy->getVectorStarts();
     3661  const int * columnLengthIn = columnCopy->getVectorLengths();
     3662  const double * elementByColumnIn = columnCopy->getElements();
     3663  int * columnLength = new int [numberColumns_];
     3664  // May just be that is a network - worth checking
     3665  bool isNetworkAlready=true;
     3666  bool trueNetwork=true;
     3667  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3668    double product =1.0;
     3669    int n=0;
     3670    for (CoinBigIndex j=columnStartIn[iColumn];j<columnStartIn[iColumn]+columnLengthIn[iColumn];j++) {
     3671      iRow = mapping[rowIn[j]];
     3672      if (iRow>=0) {
     3673        n++;
     3674        product *= elementByColumnIn[j];
     3675      }
     3676    }
     3677    if (n>=2) {
     3678      if (product!=-1.0||n>2)
     3679        isNetworkAlready=false;
     3680    } else if (n==1) {
     3681      trueNetwork=false;
     3682    }
     3683    columnLength[iColumn]=n;
     3684  }
     3685  if (!isNetworkAlready) {
     3686    // For sorting
     3687    double * count = new double [numberRows_];
     3688    int * which = new int [numberRows_];
     3689    int numberLast=-1;
     3690      // Count for columns
     3691    char * columnCount = new char[numberColumns_];
     3692    memset(columnCount,0,numberColumns_);
     3693    char * currentColumnCount = new char[numberColumns_];
     3694    // Now do main loop
     3695    while (numberIn>numberLast) {
     3696      numberLast=numberIn;
     3697      int numberLeft = 0;
     3698      for (iRow=0;iRow<numberEligible;iRow++) {
     3699        if (rotate[iRow]==0) {
     3700          which[numberLeft]=iRow;
     3701          int merit=0;
     3702          bool OK=true;
     3703          bool reflectionOK=true;
     3704          for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     3705            iColumn = column[j];
     3706            int iCount=columnCount[iColumn];
     3707            int absCount=CoinAbs(iCount);
     3708            if (absCount<2) {
     3709              merit = CoinMax(columnLength[iColumn]-absCount-1,merit);
     3710              if (elementByRow[j]==iCount)
     3711                OK=false;
     3712              else if (elementByRow[j]==-iCount)
     3713                reflectionOK=false;
     3714            } else {
     3715              merit =-2;
     3716              break;
     3717            }
     3718          }
     3719          if (merit>-2&&(OK||reflectionOK)) {
     3720            assert (!OK||!reflectionOK||!numberIn);
     3721            //if (!numberLast) merit=1;
     3722            count[numberLeft++]=(rowStart[iRow+1]-rowStart[iRow]-1)*((double)merit);
     3723            if (OK)
     3724              rotate[iRow]=0;
     3725            else
     3726              rotate[iRow]=1;
     3727          } else {
     3728            // no good
     3729            rotate[iRow]=-1;
     3730          }
     3731        }
     3732      }
     3733      CoinSort_2(count,count+numberLeft,which);
     3734      // Get G
     3735      memset(currentColumnCount,0,numberColumns_);
     3736      for (iRow=0;iRow<numberLeft;iRow++) {
     3737        int jRow = which[iRow];
     3738        bool possible=true;
     3739        for (int i=0;i<numberIn;i++) {
     3740          for (CoinBigIndex j=rowStart[jRow];j<rowStart[jRow+1];j++) {
     3741            if (currentColumnCount[column[j]]) {
     3742              possible=false;
     3743              break;
     3744            }
     3745          }
     3746        }
     3747        if (possible) {
     3748          rotate[jRow]+=2;
     3749          eligible[numberIn++]=jRow;
     3750          char multiplier = (rotate[jRow]==2) ? 1 : -1;
     3751          for (CoinBigIndex j=rowStart[jRow];j<rowStart[jRow+1];j++) {
     3752            iColumn = column[j];
     3753            currentColumnCount[iColumn]++;
     3754            int iCount=columnCount[iColumn];
     3755            int absCount=CoinAbs(iCount);
     3756            if (!absCount) {
     3757              columnCount[iColumn]=elementByRow[j]*multiplier;
     3758            } else {
     3759              columnCount[iColumn]=2;
     3760            }
     3761          }
     3762        }
     3763      }
     3764    }
     3765    for (iRow=0;iRow<numberIn;iRow++) {
     3766      int kRow = eligible[iRow];
     3767      assert (rotate[kRow]>=2);
     3768    }
     3769    trueNetwork=true;
     3770    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3771      if (CoinAbs(columnCount[iColumn])==1) {
     3772        trueNetwork=false;
     3773        break;
     3774      }
     3775    }
     3776    delete [] currentColumnCount;
     3777    delete [] columnCount;
     3778    delete [] which;
     3779    delete [] count;
     3780  } else {
     3781    numberIn=numberEligible;
     3782    for (iRow=0;iRow<numberRows_;iRow++) {
     3783      int kRow = mapping[iRow];
     3784      if (kRow>=0) {
     3785        rotate[kRow]=2;
     3786      }
     3787    }
     3788  }
     3789  if (!trueNetwork)
     3790    numberIn= - numberIn;
     3791  delete [] column;
     3792  delete [] rowStart;
     3793  delete [] elementByRow;
     3794  delete [] columnLength;
     3795  // redo rotate
     3796  char * rotate2 = CoinCopyOfArray(rotate,numberEligible);
     3797  for (iRow=0;iRow<numberRows_;iRow++) {
     3798    int kRow = mapping[iRow];
     3799    if (kRow>=0) {
     3800      int iState=rotate2[kRow];
     3801      if (iState>1)
     3802        iState -=2;
     3803      else
     3804        iState = -1;
     3805      rotate[iRow]=iState;
     3806    } else {
     3807      rotate[iRow]=-1;
     3808    }
     3809  }
     3810  delete [] rotate2;
     3811  delete [] eligible;
     3812  delete [] mapping;
     3813  return numberIn;
     3814}
    35833815//#############################################################################
    35843816// Constructors / Destructor / Assignment
  • stable/1.6/Clp/src/ClpModel.hpp

    r1060 r1120  
    276276  void setColumnName(int colIndex, std::string & name) ;
    277277#endif
     278  /** Find a network subset.
     279      rotate array should be numberRows.  On output
     280      -1 not in network
     281       0 in network as is
     282       1 in network with signs swapped
     283      Returns number of network rows
     284  */
     285  int findNetwork(char * rotate, double fractionNeeded=0.75);
    278286  /** This creates a coinModel object
    279287  */
Note: See TracChangeset for help on using the changeset viewer.