Changeset 425 for branches


Ignore:
Timestamp:
Sep 15, 2006 4:57:25 PM (13 years ago)
Author:
forrest
Message:

for sos

Location:
branches/devel/Cbc/examples
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Cbc/examples/CbcBranchLink.cpp

    r424 r425  
    1515#include "CbcBranchLink.hpp"
    1616#include "CoinError.hpp"
     17#include "CoinPackedMatrix.hpp"
    1718
    1819// Default Constructor
     
    2223    numberMembers_(0),
    2324    numberLinks_(0),
    24     which_(NULL)
     25    which_(NULL),
     26    sosType_(1)
    2527{
    2628}
     
    3234    numberMembers_(numberMembers),
    3335    numberLinks_(numberLinks),
    34     which_(NULL)
     36    which_(NULL),
     37    sosType_(1)
    3538{
    3639  id_=identifier;
     
    6164// Useful constructor (which are indices)
    6265CbcLink::CbcLink (CbcModel * model,  int numberMembers,
    63            int numberLinks, const int * which , const double * weights, int identifier)
     66           int numberLinks, int sosType, const int * which , const double * weights, int identifier)
    6467  : CbcObject(model),
    6568    numberMembers_(numberMembers),
    6669    numberLinks_(numberLinks),
    67     which_(NULL)
     70    which_(NULL),
     71    sosType_(sosType)
    6872{
    6973  id_=identifier;
     
    98102  numberMembers_ = rhs.numberMembers_;
    99103  numberLinks_ = rhs.numberLinks_;
     104  sosType_ = rhs.sosType_;
    100105  if (numberMembers_) {
    101106    weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_);
     
    124129    numberMembers_ = rhs.numberMembers_;
    125130    numberLinks_ = rhs.numberLinks_;
     131    sosType_ = rhs.sosType_;
    126132    if (numberMembers_) {
    127133      weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_);
     
    180186#endif
    181187        }
    182         weight += weights_[j]*CoinMin(value,upper[iColumn]);
     188        value = CoinMin(value,upper[iColumn]);
     189        weight += weights_[j]*value;
    183190        if (firstNonZero<0)
    184191          firstNonZero=j;
     
    188195    base += numberLinks_;
    189196  }
     197  double valueInfeasibility;
    190198  preferredWay=1;
    191   if (lastNonZero-firstNonZero>=1) {
     199  if (lastNonZero-firstNonZero>=sosType_) {
    192200    // find where to branch
    193201    assert (sum>0.0);
    194202    weight /= sum;
    195     double value = lastNonZero-firstNonZero+1;
    196     value *= 0.5/((double) numberMembers_);
    197     return value;
     203    valueInfeasibility = lastNonZero-firstNonZero+1;
     204    valueInfeasibility *= 0.5/((double) numberMembers_);
     205    //#define DISTANCE
     206#ifdef DISTANCE
     207    assert (sosType_==1); // code up
     208    /* may still be satisfied.
     209       For LOS type 2 we might wish to move coding around
     210       and keep initial info in model_ for speed
     211    */
     212    int iWhere;
     213    bool possible=false;
     214    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
     215      if (fabs(weight-weights_[iWhere])<1.0e-8) {
     216        possible=true;
     217        break;
     218      }
     219    }
     220    if (possible) {
     221      // One could move some of this (+ arrays) into model_
     222      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
     223      const double * element = matrix->getMutableElements();
     224      const int * row = matrix->getIndices();
     225      const CoinBigIndex * columnStart = matrix->getVectorStarts();
     226      const int * columnLength = matrix->getVectorLengths();
     227      const double * rowSolution = solver->getRowActivity();
     228      const double * rowLower = solver->getRowLower();
     229      const double * rowUpper = solver->getRowUpper();
     230      int numberRows = matrix->getNumRows();
     231      double * array = new double [numberRows];
     232      CoinZeroN(array,numberRows);
     233      int * which = new int [numberRows];
     234      int n=0;
     235      int base=numberLinks_*firstNonZero;
     236      for (j=firstNonZero;j<=lastNonZero;j++) {
     237        for (int k=0;k<numberLinks_;k++) {
     238          int iColumn = which_[base+k];
     239          double value = CoinMax(0.0,solution[iColumn]);
     240          if (value>integerTolerance&&upper[iColumn]) {
     241            value = CoinMin(value,upper[iColumn]);
     242            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
     243              int iRow = row[j];
     244              double a = array[iRow];
     245              if (a) {
     246                a += value*element[j];
     247                if (!a)
     248                  a = 1.0e-100;
     249              } else {
     250                which[n++]=iRow;
     251                a=value*element[j];
     252                assert (a);
     253              }
     254              array[iRow]=a;
     255            }
     256          }
     257        }
     258        base += numberLinks_;
     259      }
     260      base=numberLinks_*iWhere;
     261      for (int k=0;k<numberLinks_;k++) {
     262        int iColumn = which_[base+k];
     263        const double value = 1.0;
     264        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
     265          int iRow = row[j];
     266          double a = array[iRow];
     267          if (a) {
     268            a -= value*element[j];
     269            if (!a)
     270              a = 1.0e-100;
     271          } else {
     272            which[n++]=iRow;
     273            a=-value*element[j];
     274            assert (a);
     275          }
     276          array[iRow]=a;
     277        }
     278      }
     279      for (j=0;j<n;j++) {
     280        int iRow = which[j];
     281        // moving to point will increase row solution by this
     282        double distance = array[iRow];
     283        if (distance>1.0e-8) {
     284          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
     285            possible=false;
     286            break;
     287          }
     288        } else if (distance<-1.0e-8) {
     289          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
     290            possible=false;
     291            break;
     292          }
     293        }
     294      }
     295      for (j=0;j<n;j++)
     296        array[which[j]]=0.0;
     297      delete [] array;
     298      delete [] which;
     299      if (possible) {
     300        valueInfeasibility=0.0;
     301        printf("possible %d %d %d\n",firstNonZero,lastNonZero,iWhere);
     302      }
     303    }
     304#endif
    198305  } else {
    199     return 0.0; // satisfied
    200   }
     306    valueInfeasibility = 0.0; // satisfied
     307  }
     308  return valueInfeasibility;
    201309}
    202310
     
    231339    base += numberLinks_;
    232340  }
    233   assert (lastNonZero-firstNonZero==0) ;
     341#ifdef DISTANCE
     342  if (lastNonZero-firstNonZero>sosType_-1) {
     343    /* may still be satisfied.
     344       For LOS type 2 we might wish to move coding around
     345       and keep initial info in model_ for speed
     346    */
     347    int iWhere;
     348    bool possible=false;
     349    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
     350      if (fabs(weight-weights_[iWhere])<1.0e-8) {
     351        possible=true;
     352        break;
     353      }
     354    }
     355    if (possible) {
     356      // One could move some of this (+ arrays) into model_
     357      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
     358      const double * element = matrix->getMutableElements();
     359      const int * row = matrix->getIndices();
     360      const CoinBigIndex * columnStart = matrix->getVectorStarts();
     361      const int * columnLength = matrix->getVectorLengths();
     362      const double * rowSolution = solver->getRowActivity();
     363      const double * rowLower = solver->getRowLower();
     364      const double * rowUpper = solver->getRowUpper();
     365      int numberRows = matrix->getNumRows();
     366      double * array = new double [numberRows];
     367      CoinZeroN(array,numberRows);
     368      int * which = new int [numberRows];
     369      int n=0;
     370      int base=numberLinks_*firstNonZero;
     371      for (j=firstNonZero;j<=lastNonZero;j++) {
     372        for (int k=0;k<numberLinks_;k++) {
     373          int iColumn = which_[base+k];
     374          double value = CoinMax(0.0,solution[iColumn]);
     375          if (value>integerTolerance&&upper[iColumn]) {
     376            value = CoinMin(value,upper[iColumn]);
     377            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
     378              int iRow = row[j];
     379              double a = array[iRow];
     380              if (a) {
     381                a += value*element[j];
     382                if (!a)
     383                  a = 1.0e-100;
     384              } else {
     385                which[n++]=iRow;
     386                a=value*element[j];
     387                assert (a);
     388              }
     389              array[iRow]=a;
     390            }
     391          }
     392        }
     393        base += numberLinks_;
     394      }
     395      base=numberLinks_*iWhere;
     396      for (int k=0;k<numberLinks_;k++) {
     397        int iColumn = which_[base+k];
     398        const double value = 1.0;
     399        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
     400          int iRow = row[j];
     401          double a = array[iRow];
     402          if (a) {
     403            a -= value*element[j];
     404            if (!a)
     405              a = 1.0e-100;
     406          } else {
     407            which[n++]=iRow;
     408            a=-value*element[j];
     409            assert (a);
     410          }
     411          array[iRow]=a;
     412        }
     413      }
     414      for (j=0;j<n;j++) {
     415        int iRow = which[j];
     416        // moving to point will increase row solution by this
     417        double distance = array[iRow];
     418        if (distance>1.0e-8) {
     419          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
     420            possible=false;
     421            break;
     422          }
     423        } else if (distance<-1.0e-8) {
     424          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
     425            possible=false;
     426            break;
     427          }
     428        }
     429      }
     430      for (j=0;j<n;j++)
     431        array[which[j]]=0.0;
     432      delete [] array;
     433      delete [] which;
     434      if (possible) {
     435        printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere);
     436        firstNonZero=iWhere;
     437        lastNonZero=iWhere;
     438      }
     439    }
     440  }
     441#else
     442  assert (lastNonZero-firstNonZero<sosType_) ;
     443#endif
    234444  base=0;
    235445  for (j=0;j<firstNonZero;j++) {
     
    288498    base += numberLinks_;
    289499  }
    290   assert (lastNonZero-firstNonZero>=1) ;
     500  assert (lastNonZero-firstNonZero>=sosType_) ;
    291501  // find where to branch
    292502  assert (sum>0.0);
     
    297507    if (weight<weights_[iWhere+1])
    298508      break;
    299   separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
     509  if (sosType_==1) {
     510    // SOS 1
     511    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
     512  } else {
     513    // SOS 2
     514    if (iWhere==firstNonFixed)
     515      iWhere++;;
     516    if (iWhere==lastNonFixed-1)
     517      iWhere = lastNonFixed-2;
     518    separator = weights_[iWhere+1];
     519  }
    300520  // create object
    301521  CbcBranchingObject * branch;
  • branches/devel/Cbc/examples/CbcBranchLink.hpp

    r424 r425  
    3434  */
    3535  CbcLink (CbcModel * model, int numberMembers,
    36            int numberLinks, const int * which,
     36           int numberLinks, int typeSOS, const int * which,
    3737           const double * weights, int setNumber);
    3838 
     
    8585  /// Members
    8686  int * which_;
     87  /// Type 1 or 2
     88  int sosType_;
    8789};
    8890/** Branching object for Special ordered sets
Note: See TracChangeset for help on using the changeset viewer.