Changeset 912


Ignore:
Timestamp:
Apr 10, 2008 1:55:10 PM (11 years ago)
Author:
ladanyi
Message:

Incorporated changes from branches/heur

Location:
trunk/Cbc/src
Files:
22 edited
2 copied

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcBranchActual.cpp

    r904 r912  
    1111//#define CBC_DEBUG
    1212
     13#include "CoinTypes.hpp"
    1314#include "OsiSolverInterface.hpp"
    1415#include "OsiSolverBranch.hpp"
     
    1819#include "CoinSort.hpp"
    1920#include "CoinError.hpp"
     21
     22//##############################################################################
    2023
    2124// Default Constructor
     
    332335  return branch;
    333336}
     337
     338//##############################################################################
    334339
    335340// Default Constructor
     
    709714}
    710715
     716//##############################################################################
     717
    711718/** Default Constructor
    712719
     
    939946  return NULL;
    940947}
     948
     949//##############################################################################
     950
    941951// Default Constructor
    942952CbcIntegerBranchingObject::CbcIntegerBranchingObject()
     
    13321342}
    13331343
     1344/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     1345    same type and must have the same original object, but they may have
     1346    different feasible regions.
     1347    Return the appropriate CbcRangeCompare value (first argument being the
     1348    sub/superset if that's the case). In case of overlap (and if \c
     1349    replaceIfOverlap is true) replace the current branching object with one
     1350    whose feasible region is the overlap.
     1351   */
     1352CbcRangeCompare
     1353CbcIntegerBranchingObject::compareBranchingObject
     1354(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     1355{
     1356  const CbcIntegerBranchingObject* br =
     1357    dynamic_cast<const CbcIntegerBranchingObject*>(brObj);
     1358  assert(br);
     1359  double* thisBd = way_ < 0 ? down_ : up_;
     1360  const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
     1361  return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
     1362}
     1363
     1364//##############################################################################
    13341365
    13351366/** Default Constructor
     
    15691600}
    15701601
     1602//##############################################################################
     1603
    15711604// Default Constructor
    15721605CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
     
    16401673}
    16411674
     1675/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     1676    same type and must have the same original object, but they may have
     1677    different feasible regions.
     1678    Return the appropriate CbcRangeCompare value (first argument being the
     1679    sub/superset if that's the case). In case of overlap (and if \c
     1680    replaceIfOverlap is true) replace the current branching object with one
     1681    whose feasible region is the overlap.
     1682*/
     1683CbcRangeCompare
     1684CbcIntegerPseudoCostBranchingObject::compareBranchingObject
     1685(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     1686{
     1687  const CbcIntegerPseudoCostBranchingObject* br =
     1688    dynamic_cast<const CbcIntegerPseudoCostBranchingObject*>(brObj);
     1689  assert(br);
     1690  double* thisBd = way_ < 0 ? down_ : up_;
     1691  const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
     1692  return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
     1693}
     1694
     1695
     1696//##############################################################################
    16421697
    16431698// Default Constructor
     
    18141869  printf("\n");
    18151870}
    1816  
     1871
     1872static inline int
     1873CbcCompareCliques(const CbcClique* cl0, const CbcClique* cl1)
     1874{
     1875  if (cl0->cliqueType() < cl1->cliqueType()) {
     1876    return -1;
     1877  }
     1878  if (cl0->cliqueType() > cl1->cliqueType()) {
     1879    return 1;
     1880  }
     1881  if (cl0->numberMembers() != cl1->numberMembers()) {
     1882    return cl0->numberMembers() - cl1->numberMembers();
     1883  }
     1884  if (cl0->numberNonSOSMembers() != cl1->numberNonSOSMembers()) {
     1885    return cl0->numberNonSOSMembers() - cl1->numberNonSOSMembers();
     1886  }
     1887  return memcmp(cl0->members(), cl1->members(),
     1888                cl0->numberMembers() * sizeof(int));
     1889}
     1890
     1891/** Compare the original object of \c this with the original object of \c
     1892    brObj. Assumes that there is an ordering of the original objects.
     1893    This method should be invoked only if \c this and brObj are of the same
     1894    type.
     1895    Return negative/0/positive depending on whether \c this is
     1896    smaller/same/larger than the argument.
     1897*/
     1898int
     1899CbcCliqueBranchingObject::compareOriginalObject
     1900(const CbcBranchingObject* brObj) const
     1901{
     1902  const CbcCliqueBranchingObject* br =
     1903    dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
     1904  assert(br);
     1905  return CbcCompareCliques(clique_, br->clique_);
     1906}
     1907
     1908/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     1909    same type and must have the same original object, but they may have
     1910    different feasible regions.
     1911    Return the appropriate CbcRangeCompare value (first argument being the
     1912    sub/superset if that's the case). In case of overlap (and if \c
     1913    replaceIfOverlap is true) replace the current branching object with one
     1914    whose feasible region is the overlap.
     1915*/
     1916CbcRangeCompare
     1917CbcCliqueBranchingObject::compareBranchingObject
     1918(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     1919{
     1920  const CbcCliqueBranchingObject* br =
     1921    dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
     1922  assert(br);
     1923  unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
     1924  const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
     1925  const CoinUInt64 cl0 =
     1926    (static_cast<CoinUInt64>(thisMask[0]) << 32) | thisMask[1];
     1927  const CoinUInt64 cl1 =
     1928    (static_cast<CoinUInt64>(otherMask[0]) << 32) | otherMask[1];
     1929  if (cl0 == cl1) {
     1930    return CbcRangeSame;
     1931  }
     1932  const CoinUInt64 cl_intersection = (cl0 & cl1);
     1933  if (cl_intersection == cl0) {
     1934    return CbcRangeSuperset;
     1935  }
     1936  if (cl_intersection == cl1) {
     1937    return CbcRangeSubset;
     1938  }
     1939  const CoinUInt64 cl_xor = (cl0 ^ cl1);
     1940  if (cl_intersection == 0 && cl_xor == 0) {
     1941    return CbcRangeDisjoint;
     1942  }
     1943  const CoinUInt64 cl_union = (cl0 | cl1);
     1944  thisMask[0] = static_cast<unsigned int>(cl_union >> 32);
     1945  thisMask[1] = static_cast<unsigned int>(cl_union & 0xffffffff);
     1946  return CbcRangeOverlap;
     1947}
     1948
     1949//##############################################################################
     1950
    18171951// Default Constructor
    18181952CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject()
     
    20052139  printf("\n");
    20062140}
     2141
     2142/** Compare the original object of \c this with the original object of \c
     2143    brObj. Assumes that there is an ordering of the original objects.
     2144    This method should be invoked only if \c this and brObj are of the same
     2145    type.
     2146    Return negative/0/positive depending on whether \c this is
     2147    smaller/same/larger than the argument.
     2148*/
     2149int
     2150CbcLongCliqueBranchingObject::compareOriginalObject
     2151(const CbcBranchingObject* brObj) const
     2152{
     2153  const CbcLongCliqueBranchingObject* br =
     2154    dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
     2155  assert(br);
     2156  return CbcCompareCliques(clique_, br->clique_);
     2157}
     2158
     2159/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     2160    same type and must have the same original object, but they may have
     2161    different feasible regions.
     2162    Return the appropriate CbcRangeCompare value (first argument being the
     2163    sub/superset if that's the case). In case of overlap (and if \c
     2164    replaceIfOverlap is true) replace the current branching object with one
     2165    whose feasible region is the overlap.
     2166*/
     2167CbcRangeCompare
     2168CbcLongCliqueBranchingObject::compareBranchingObject
     2169(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     2170{
     2171  const CbcLongCliqueBranchingObject* br =
     2172    dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
     2173  assert(br);
     2174  const int numberMembers = clique_->numberMembers();
     2175  const int numberWords=(numberMembers+31)>>5;
     2176  unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
     2177  const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
     2178
     2179  if (memcmp(thisMask, otherMask, numberWords * sizeof(unsigned int)) == 0) {
     2180    return CbcRangeSame;
     2181  }
     2182  bool canBeSuperset = true;
     2183  bool canBeSubset = true;
     2184  int i;
     2185  for (i = numberWords-1; i >= 0 && (canBeSuperset || canBeSubset); --i) {
     2186    const unsigned int both = (thisMask[i] & otherMask[i]);
     2187    canBeSuperset &= (both == thisMask[i]);
     2188    canBeSubset &= (both == otherMask[i]);
     2189  }
     2190  if (canBeSuperset) {
     2191    return CbcRangeSuperset;
     2192  }
     2193  if (canBeSubset) {
     2194    return CbcRangeSubset;
     2195  }
     2196
     2197  for (i = numberWords-1; i >= 0; --i) {
     2198    if ((thisMask[i] ^ otherMask[i]) != 0) {
     2199      break;
     2200    }
     2201  }
     2202  if (i == -1) { // complement
     2203    return CbcRangeDisjoint;
     2204  }
     2205  // must be overlap
     2206  for (i = numberWords-1; i >= 0; --i) {
     2207    thisMask[i] |= otherMask[i];
     2208  }
     2209  return CbcRangeOverlap;
     2210}
     2211
     2212//##############################################################################
     2213
    20072214// Default Constructor
    20082215CbcSOSBranchingObject::CbcSOSBranchingObject()
    2009   :CbcBranchingObject()
     2216  :CbcBranchingObject(),
     2217   firstNonzero_(-1),
     2218   lastNonzero_(-1)
    20102219{
    20112220  set_ = NULL;
     
    20222231  set_ = set;
    20232232  separator_ = separator;
     2233  computeNonzeroRange();
    20242234}
    20252235
    20262236// Copy constructor
    2027 CbcSOSBranchingObject::CbcSOSBranchingObject ( const CbcSOSBranchingObject & rhs) :CbcBranchingObject(rhs)
     2237CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
     2238 :CbcBranchingObject(rhs),
     2239  firstNonzero_(rhs.firstNonzero_),
     2240  lastNonzero_(rhs.lastNonzero_)
    20282241{
    20292242  set_=rhs.set_;
     
    20392252    set_=rhs.set_;
    20402253    separator_ = rhs.separator_;
     2254    firstNonzero_ = rhs.firstNonzero_;
     2255    lastNonzero_ = rhs.lastNonzero_;
    20412256  }
    20422257  return *this;
     
    20532268{
    20542269}
     2270
     2271void
     2272CbcSOSBranchingObject::computeNonzeroRange()
     2273{
     2274  const int numberMembers = set_->numberMembers();
     2275  const double * weights = set_->weights();
     2276  int i = 0;
     2277  if (way_ < 0) {
     2278    for ( i=0;i<numberMembers;i++) {
     2279      if (weights[i] > separator_)
     2280        break;
     2281    }
     2282    assert (i<numberMembers);
     2283    firstNonzero_ = 0;
     2284    lastNonzero_ = i;
     2285  } else {
     2286    for ( i=0;i<numberMembers;i++) {
     2287      if (weights[i] >= separator_)
     2288        break;
     2289    }
     2290    assert (i<numberMembers);
     2291    firstNonzero_ = i;
     2292    lastNonzero_ = numberMembers;
     2293  }
     2294}
     2295
    20552296double
    20562297CbcSOSBranchingObject::branch()
     
    20852326    way_=-1;      // Swap direction
    20862327  }
     2328  computeNonzeroRange();
    20872329  return 0.0;
    20882330}
     
    21452387}
    21462388 
     2389/** Compare the original object of \c this with the original object of \c
     2390    brObj. Assumes that there is an ordering of the original objects.
     2391    This method should be invoked only if \c this and brObj are of the same
     2392    type.
     2393    Return negative/0/positive depending on whether \c this is
     2394    smaller/same/larger than the argument.
     2395*/
     2396int
     2397CbcSOSBranchingObject::compareOriginalObject
     2398(const CbcBranchingObject* brObj) const
     2399{
     2400  const CbcSOSBranchingObject* br =
     2401    dynamic_cast<const CbcSOSBranchingObject*>(brObj);
     2402  assert(br);
     2403  const CbcSOS* s0 = set_;
     2404  const CbcSOS* s1 = br->set_;
     2405  if (s0->sosType() != s1->sosType()) {
     2406    return s0->sosType() - s1->sosType();
     2407  }
     2408  if (s0->numberMembers() != s1->numberMembers()) {
     2409    return s0->numberMembers() - s1->numberMembers();
     2410  }
     2411  const int memberCmp = memcmp(s0->members(), s1->members(),
     2412                               s0->numberMembers() * sizeof(int));
     2413  if (memberCmp != 0) {
     2414    return memberCmp;
     2415  }
     2416  return memcmp(s0->weights(), s1->weights(),
     2417                s0->numberMembers() * sizeof(double));
     2418}
     2419
     2420/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     2421    same type and must have the same original object, but they may have
     2422    different feasible regions.
     2423    Return the appropriate CbcRangeCompare value (first argument being the
     2424    sub/superset if that's the case). In case of overlap (and if \c
     2425    replaceIfOverlap is true) replace the current branching object with one
     2426    whose feasible region is the overlap.
     2427*/
     2428CbcRangeCompare
     2429CbcSOSBranchingObject::compareBranchingObject
     2430(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     2431{
     2432  const CbcSOSBranchingObject* br =
     2433    dynamic_cast<const CbcSOSBranchingObject*>(brObj);
     2434  assert(br);
     2435  if (firstNonzero_ < br->firstNonzero_) {
     2436    if (lastNonzero_ >= br->lastNonzero_) {
     2437      return CbcRangeSuperset;
     2438    } else if (lastNonzero_ <= br->firstNonzero_) {
     2439      return CbcRangeDisjoint;
     2440    } else {
     2441      // overlap
     2442      if (replaceIfOverlap) {
     2443        firstNonzero_ = br->firstNonzero_;
     2444      }
     2445      return CbcRangeOverlap;
     2446    }
     2447  } else if (firstNonzero_ > br->firstNonzero_) {
     2448    if (lastNonzero_ <= br->lastNonzero_) {
     2449      return CbcRangeSubset;
     2450    } else if (firstNonzero_ >= br->lastNonzero_) {
     2451      return CbcRangeDisjoint;
     2452    } else {
     2453      // overlap
     2454      if (replaceIfOverlap) {
     2455        lastNonzero_ = br->lastNonzero_;
     2456      }
     2457      return CbcRangeOverlap;
     2458    }
     2459  } else {
     2460    if (lastNonzero_ == br->lastNonzero_) {
     2461      return CbcRangeSame;
     2462    }
     2463    return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
     2464  }
     2465  return CbcRangeSame; // fake return
     2466}
     2467
     2468//##############################################################################
     2469
    21472470// Default Constructor
    21482471CbcBranchDefaultDecision::CbcBranchDefaultDecision()
     
    25552878}
    25562879
     2880//##############################################################################
     2881
    25572882// Default Constructor
    25582883CbcFollowOn::CbcFollowOn ()
     
    28933218  return branch;
    28943219}
     3220
     3221//##############################################################################
     3222
    28953223// Default Constructor
    28963224CbcFixingBranchingObject::CbcFixingBranchingObject()
     
    30133341  printf("\n");
    30143342}
     3343
     3344/** Compare the original object of \c this with the original object of \c
     3345    brObj. Assumes that there is an ordering of the original objects.
     3346    This method should be invoked only if \c this and brObj are of the same
     3347    type.
     3348    Return negative/0/positive depending on whether \c this is
     3349    smaller/same/larger than the argument.
     3350*/
     3351int
     3352CbcFixingBranchingObject::compareOriginalObject
     3353(const CbcBranchingObject* brObj) const
     3354{
     3355  throw("must implement");
     3356}
     3357
     3358/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     3359    same type and must have the same original object, but they may have
     3360    different feasible regions.
     3361    Return the appropriate CbcRangeCompare value (first argument being the
     3362    sub/superset if that's the case). In case of overlap (and if \c
     3363    replaceIfOverlap is true) replace the current branching object with one
     3364    whose feasible region is the overlap.
     3365   */
     3366CbcRangeCompare
     3367CbcFixingBranchingObject::compareBranchingObject
     3368(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     3369{
     3370  const CbcFixingBranchingObject* br =
     3371    dynamic_cast<const CbcFixingBranchingObject*>(brObj);
     3372  assert(br);
     3373  // If two FixingBranchingObject's have the same base object then it's pretty
     3374  // much guaranteed
     3375  throw("must implement");
     3376}
     3377
     3378//##############################################################################
     3379
    30153380// Default Constructor
    30163381CbcNWay::CbcNWay ()
     
    32623627}
    32633628 
     3629//##############################################################################
     3630
    32643631// Default Constructor
    32653632CbcNWayBranchingObject::CbcNWayBranchingObject()
     
    33833750  printf("\n");
    33843751}
     3752
     3753/** Compare the original object of \c this with the original object of \c
     3754    brObj. Assumes that there is an ordering of the original objects.
     3755    This method should be invoked only if \c this and brObj are of the same
     3756    type.
     3757    Return negative/0/positive depending on whether \c this is
     3758    smaller/same/larger than the argument.
     3759*/
     3760int
     3761CbcNWayBranchingObject::compareOriginalObject
     3762(const CbcBranchingObject* brObj) const
     3763{
     3764  throw("must implement");
     3765}
     3766
     3767/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     3768    same type and must have the same original object, but they may have
     3769    different feasible regions.
     3770    Return the appropriate CbcRangeCompare value (first argument being the
     3771    sub/superset if that's the case). In case of overlap (and if \c
     3772    replaceIfOverlap is true) replace the current branching object with one
     3773    whose feasible region is the overlap.
     3774*/
     3775CbcRangeCompare
     3776CbcNWayBranchingObject::compareBranchingObject
     3777(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     3778{
     3779  throw("must implement");
     3780}
     3781
     3782//##############################################################################
    33853783
    33863784// Default Constructor
     
    35503948}
    35513949
     3950//##############################################################################
     3951
    35523952// Default Constructor
    35533953CbcDummyBranchingObject::CbcDummyBranchingObject(CbcModel * model)
     
    35993999  printf("Dummy branch\n");
    36004000}
     4001
     4002/** Compare the original object of \c this with the original object of \c
     4003    brObj. Assumes that there is an ordering of the original objects.
     4004    This method should be invoked only if \c this and brObj are of the same
     4005    type.
     4006    Return negative/0/positive depending on whether \c this is
     4007    smaller/same/larger than the argument.
     4008*/
     4009int
     4010CbcDummyBranchingObject::compareOriginalObject
     4011(const CbcBranchingObject* brObj) const
     4012{
     4013  throw("must implement");
     4014}
     4015
     4016/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     4017    same type and must have the same original object, but they may have
     4018    different feasible regions.
     4019    Return the appropriate CbcRangeCompare value (first argument being the
     4020    sub/superset if that's the case). In case of overlap (and if \c
     4021    replaceIfOverlap is true) replace the current branching object with one
     4022    whose feasible region is the overlap.
     4023*/
     4024CbcRangeCompare
     4025CbcDummyBranchingObject::compareBranchingObject
     4026(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     4027{
     4028  throw("must implement");
     4029}
     4030
  • trunk/Cbc/src/CbcBranchActual.hpp

    r848 r912  
    449449  */
    450450  virtual double branch();
     451
     452#if 0
     453  // No need to override. Default works fine.
     454  /** Reset every information so that the branching object appears to point to
     455      the previous child. This method does not need to modify anything in any
     456      solver. */
     457  virtual void previousBranch();
     458#endif
    451459
    452460  using CbcBranchingObject::print ;
     
    487495#endif
    488496
     497  /** Return the type (an integer identifier) of \c this */
     498  virtual int type() const { return 100; }
     499
     500  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     501      same type and must have the same original object, but they may have
     502      different feasible regions.
     503      Return the appropriate CbcRangeCompare value (first argument being the
     504      sub/superset if that's the case). In case of overlap (and if \c
     505      replaceIfOverlap is true) replace the current branching object with one
     506      whose feasible region is the overlap.
     507   */
     508  virtual CbcRangeCompare compareBranchingObject
     509  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     510
    489511protected:
    490512  /// Lower [0] and upper [1] bounds for the down arm (way_ = -1)
     
    661683  inline void setChangeInGuessed(double value)
    662684  { changeInGuessed_=value;}
     685
     686  /** Return the type (an integer identifier) of \c this */
     687  virtual int type() const { return 101; }
     688
     689  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     690      same type and must have the same original object, but they may have
     691      different feasible regions.
     692      Return the appropriate CbcRangeCompare value (first argument being the
     693      sub/superset if that's the case). In case of overlap (and if \c
     694      replaceIfOverlap is true) replace the current branching object with one
     695      whose feasible region is the overlap.
     696   */
     697  virtual CbcRangeCompare compareBranchingObject
     698  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     699
    663700protected:
    664701  /// Change in guessed objective value for next branch
     
    704741  virtual double branch();
    705742
     743#if 0
     744  // No need to override. Default works fine.
     745  /** Reset every information so that the branching object appears to point to
     746      the previous child. This method does not need to modify anything in any
     747      solver. */
     748  virtual void previousBranch();
     749#endif
     750
    706751  using CbcBranchingObject::print ;
    707752  /** \brief Print something about branch - only if log level high
    708753  */
    709754  virtual void print();
     755
     756  /** Return the type (an integer identifier) of \c this */
     757  virtual int type() const { return 102; }
     758
     759  /** Compare the original object of \c this with the original object of \c
     760      brObj. Assumes that there is an ordering of the original objects.
     761      This method should be invoked only if \c this and brObj are of the same
     762      type.
     763      Return negative/0/positive depending on whether \c this is
     764      smaller/same/larger than the argument.
     765  */
     766  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
     767
     768  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     769      same type and must have the same original object, but they may have
     770      different feasible regions.
     771      Return the appropriate CbcRangeCompare value (first argument being the
     772      sub/superset if that's the case). In case of overlap (and if \c
     773      replaceIfOverlap is true) replace the current branching object with one
     774      whose feasible region is the overlap.
     775   */
     776  virtual CbcRangeCompare compareBranchingObject
     777  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     778
    710779private:
    711780  /// data
     
    751820  virtual double branch();
    752821
     822#if 0
     823  // No need to override. Default works fine.
     824  /** Reset every information so that the branching object appears to point to
     825      the previous child. This method does not need to modify anything in any
     826      solver. */
     827  virtual void previousBranch();
     828#endif
     829
    753830  using CbcBranchingObject::print ;
    754831  /** \brief Print something about branch - only if log level high
    755832  */
    756833  virtual void print();
     834
     835  /** Return the type (an integer identifier) of \c this */
     836  virtual int type() const { return 103; }
     837
     838  /** Compare the original object of \c this with the original object of \c
     839      brObj. Assumes that there is an ordering of the original objects.
     840      This method should be invoked only if \c this and brObj are of the same
     841      type.
     842      Return negative/0/positive depending on whether \c this is
     843      smaller/same/larger than the argument.
     844  */
     845  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
     846
     847  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     848      same type and must have the same original object, but they may have
     849      different feasible regions.
     850      Return the appropriate CbcRangeCompare value (first argument being the
     851      sub/superset if that's the case). In case of overlap (and if \c
     852      replaceIfOverlap is true) replace the current branching object with one
     853      whose feasible region is the overlap.
     854   */
     855  virtual CbcRangeCompare compareBranchingObject
     856  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     857
    757858private:
    758859  /// data
     
    797898  virtual double branch();
    798899
     900  /** Reset every information so that the branching object appears to point to
     901      the previous child. This method does not need to modify anything in any
     902      solver. */
     903  virtual void previousBranch() {
     904    CbcBranchingObject::previousBranch();
     905    computeNonzeroRange();
     906  }
     907
    799908  using CbcBranchingObject::print ;
    800909  /** \brief Print something about branch - only if log level high
    801910  */
    802911  virtual void print();
     912
     913  /** Return the type (an integer identifier) of \c this */
     914  virtual int type() const { return 104; }
     915
     916  /** Compare the original object of \c this with the original object of \c
     917      brObj. Assumes that there is an ordering of the original objects.
     918      This method should be invoked only if \c this and brObj are of the same
     919      type.
     920      Return negative/0/positive depending on whether \c this is
     921      smaller/same/larger than the argument.
     922  */
     923  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
     924
     925  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     926      same type and must have the same original object, but they may have
     927      different feasible regions.
     928      Return the appropriate CbcRangeCompare value (first argument being the
     929      sub/superset if that's the case). In case of overlap (and if \c
     930      replaceIfOverlap is true) replace the current branching object with one
     931      whose feasible region is the overlap.
     932   */
     933  virtual CbcRangeCompare compareBranchingObject
     934  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     935
     936  /** Fill out the \c firstNonzero_ and \c lastNonzero_ data members */
     937  void computeNonzeroRange();
     938
    803939private:
    804940  /// data
     
    806942  /// separator
    807943  double separator_;
     944  /** The following two members describe the range in the members_ of the
     945      original object that whose upper bound is not fixed to 0. This is not
     946      necessary for Cbc to function correctly, this is there for heuristics so
     947      that separate branching decisions on the same object can be pooled into
     948      one branching object. */
     949  int firstNonzero_;
     950  int lastNonzero_;
    808951};
    809952
     
    840983  /// Does next branch and updates state
    841984  virtual double branch();
     985
     986#if 0
     987  // FIXME: what do we need to do here?
     988  /** Reset every information so that the branching object appears to point to
     989      the previous child. This method does not need to modify anything in any
     990      solver. */
     991  virtual void previousBranch();
     992#endif
    842993
    843994  using CbcBranchingObject::print ;
     
    8521003  virtual bool twoWay() const
    8531004  { return false;}
     1005
     1006  /** Return the type (an integer identifier) of \c this */
     1007  virtual int type() const { return 105; }
     1008
     1009  /** Compare the original object of \c this with the original object of \c
     1010      brObj. Assumes that there is an ordering of the original objects.
     1011      This method should be invoked only if \c this and brObj are of the same
     1012      type.
     1013      Return negative/0/positive depending on whether \c this is
     1014      smaller/same/larger than the argument.
     1015  */
     1016  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
     1017
     1018  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     1019      same type and must have the same original object, but they may have
     1020      different feasible regions.
     1021      Return the appropriate CbcRangeCompare value (first argument being the
     1022      sub/superset if that's the case). In case of overlap (and if \c
     1023      replaceIfOverlap is true) replace the current branching object with one
     1024      whose feasible region is the overlap.
     1025   */
     1026  virtual CbcRangeCompare compareBranchingObject
     1027  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     1028
    8541029private:
    8551030  /// order of branching - points back to CbcNWay
     
    10371212  virtual double branch();
    10381213
     1214#if 0
     1215  // No need to override. Default works fine.
     1216  /** Reset every information so that the branching object appears to point to
     1217      the previous child. This method does not need to modify anything in any
     1218      solver. */
     1219  virtual void previousBranch();
     1220#endif
     1221
    10391222  using CbcBranchingObject::print ;
    10401223  /** \brief Print something about branch - only if log level high
    10411224  */
    10421225  virtual void print();
     1226
     1227  /** Return the type (an integer identifier) of \c this */
     1228  virtual int type() const { return 106; }
     1229
     1230  /** Compare the original object of \c this with the original object of \c
     1231      brObj. Assumes that there is an ordering of the original objects.
     1232      This method should be invoked only if \c this and brObj are of the same
     1233      type.
     1234      Return negative/0/positive depending on whether \c this is
     1235      smaller/same/larger than the argument.
     1236  */
     1237  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
     1238
     1239  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     1240      same type and must have the same original object, but they may have
     1241      different feasible regions.
     1242      Return the appropriate CbcRangeCompare value (first argument being the
     1243      sub/superset if that's the case). In case of overlap (and if \c
     1244      replaceIfOverlap is true) replace the current branching object with one
     1245      whose feasible region is the overlap.
     1246   */
     1247  virtual CbcRangeCompare compareBranchingObject
     1248  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     1249
    10431250private:
    10441251  /// data
     
    11351342  virtual double branch();
    11361343
     1344#if 0
     1345  // No need to override. Default works fine.
     1346  /** Reset every information so that the branching object appears to point to
     1347      the previous child. This method does not need to modify anything in any
     1348      solver. */
     1349  virtual void previousBranch();
     1350#endif
     1351
    11371352  using CbcBranchingObject::print ;
    11381353  /** \brief Print something about branch - only if log level high
     
    11401355  virtual void print();
    11411356
     1357  /** Return the type (an integer identifier) of \c this */
     1358  virtual int type() const { return 107; }
     1359
     1360  /** Compare the original object of \c this with the original object of \c
     1361      brObj. Assumes that there is an ordering of the original objects.
     1362      This method should be invoked only if \c this and brObj are of the same
     1363      type.
     1364      Return negative/0/positive depending on whether \c this is
     1365      smaller/same/larger than the argument.
     1366  */
     1367  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
     1368
     1369  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     1370      same type and must have the same original object, but they may have
     1371      different feasible regions.
     1372      Return the appropriate CbcRangeCompare value (first argument being the
     1373      sub/superset if that's the case). In case of overlap (and if \c
     1374      replaceIfOverlap is true) replace the current branching object with one
     1375      whose feasible region is the overlap.
     1376   */
     1377  virtual CbcRangeCompare compareBranchingObject
     1378  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
     1379
    11421380};
    11431381
  • trunk/Cbc/src/CbcBranchBase.hpp

    r833 r912  
    1616class OsiChooseVariable;
    1717class CbcObjectUpdateData;
     18
     19//#############################################################################
     20
     21enum CbcRangeCompare {
     22  CbcRangeSame,
     23  CbcRangeDisjoint,
     24  CbcRangeSubset,
     25  CbcRangeSuperset,
     26  CbcRangeOverlap
     27};
    1828
    1929//#############################################################################
     
    299309  { return branch();}
    300310
     311  /** Reset every information so that the branching object appears to point to
     312      the previous child. This method does not need to modify anything in any
     313      solver. */
     314  virtual void previousBranch() {
     315    assert(branchIndex_ > 0);
     316    branchIndex_--;
     317    way_ = -way_;
     318  }
     319
    301320  using OsiBranchingObject::print ;
    302321  /** \brief Print something about branch - only if log level high
     
    348367  inline void setOriginalObject(CbcObject * object)
    349368  {originalCbcObject_=object;}
     369
     370  // Methods used in heuristics
     371 
     372  /** Return the type (an integer identifier) of \c this */
     373  virtual int type() const = 0;
     374
     375  /** Compare the original object of \c this with the original object of \c
     376      brObj. Assumes that there is an ordering of the original objects.
     377      This method should be invoked only if \c this and brObj are of the same
     378      type.
     379      Return negative/0/positive depending on whether \c this is
     380      smaller/same/larger than the argument.
     381  */
     382  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const
     383  {
     384    const CbcBranchingObject* br=dynamic_cast<const CbcBranchingObject*>(brObj);
     385    return variable() - br->variable();
     386  }
     387
     388  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     389      same type and must have the same original object, but they may have
     390      different feasible regions.
     391      Return the appropriate CbcRangeCompare value (first argument being the
     392      sub/superset if that's the case). In case of overlap (and if \c
     393      replaceIfOverlap is true) replace the current branching object with one
     394      whose feasible region is the overlap.
     395   */
     396  virtual CbcRangeCompare compareBranchingObject
     397  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false) = 0;
    350398
    351399protected:
     
    555603};
    556604
     605//##############################################################################
     606
     607/** Compare two ranges. The two bounds arrays are both of size two and
     608    describe closed intervals. Return the appropriate CbcRangeCompare value
     609    (first argument being the sub/superset if that's the case). In case of
     610    overlap (and if \c replaceIfOverlap is true) replace the content of thisBd
     611    with the intersection of the ranges.
     612*/
     613static inline CbcRangeCompare
     614CbcCompareRanges(double* thisBd, const double* otherBd,
     615                 const bool replaceIfOverlap)
     616{
     617  const double lbDiff = thisBd[0] - otherBd[0];
     618  if (lbDiff < 0) { // lb of this < lb of other
     619    if (thisBd[1] >= otherBd[1]) { // ub of this >= ub of other
     620      return CbcRangeSuperset;
     621    } else if (thisBd[1] < otherBd[0]) {
     622      return CbcRangeDisjoint;
     623    } else {
     624      // overlap
     625      if (replaceIfOverlap) {
     626        thisBd[0] = otherBd[0];
     627      }
     628      return CbcRangeOverlap;
     629    }
     630  } else if (lbDiff > 0) { // lb of this > lb of other
     631    if (thisBd[1] <= otherBd[1]) { // ub of this <= ub of other
     632      return CbcRangeSubset;
     633    } else if (thisBd[0] > otherBd[1]) {
     634      return CbcRangeDisjoint;
     635    } else {
     636      // overlap
     637      if (replaceIfOverlap) {
     638        thisBd[1] = otherBd[1];
     639      }
     640      return CbcRangeOverlap;
     641    }
     642  } else { // lb of this == lb of other
     643    if (thisBd[1] == otherBd[1]) {
     644      return CbcRangeSame;
     645    }
     646    return thisBd[1] < otherBd[1] ? CbcRangeSubset : CbcRangeSuperset;
     647  }
     648
     649  return CbcRangeSame; // fake return
     650
     651}
     652
     653//#############################################################################
     654
    557655#endif
  • trunk/Cbc/src/CbcBranchCut.cpp

    r904 r912  
    285285  return false;
    286286}
     287
     288/** Compare the original object of \c this with the original object of \c
     289    brObj. Assumes that there is an ordering of the original objects.
     290    This method should be invoked only if \c this and brObj are of the same
     291    type.
     292    Return negative/0/positive depending on whether \c this is
     293    smaller/same/larger than the argument.
     294*/
     295int
     296CbcCutBranchingObject::compareOriginalObject
     297(const CbcBranchingObject* brObj) const
     298{
     299  const CbcCutBranchingObject* br =
     300    dynamic_cast<const CbcCutBranchingObject*>(brObj);
     301  assert(br);
     302  const OsiRowCut& r0 = way_ == -1 ? down_ : up_;
     303  const OsiRowCut& r1 = br->way_ == -1 ? br->down_ : br->up_;
     304  return r0.row().compare(r1.row());
     305}
     306
     307/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     308    same type and must have the same original object, but they may have
     309    different feasible regions.
     310    Return the appropriate CbcRangeCompare value (first argument being the
     311    sub/superset if that's the case). In case of overlap (and if \c
     312    replaceIfOverlap is true) replace the current branching object with one
     313    whose feasible region is the overlap.
     314*/
     315
     316CbcRangeCompare
     317CbcCutBranchingObject::compareBranchingObject
     318(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     319{
     320  const CbcCutBranchingObject* br =
     321    dynamic_cast<const CbcCutBranchingObject*>(brObj);
     322  assert(br);
     323  OsiRowCut& r0 = way_ == -1 ? down_ : up_;
     324  const OsiRowCut& r1 = br->way_ == -1 ? br->down_ : br->up_;
     325  double thisBd[2];
     326  thisBd[0] = r0.lb();
     327  thisBd[1] = r0.ub();
     328  double otherBd[2];
     329  otherBd[0] = r1.lb();
     330  otherBd[1] = r1.ub();
     331  CbcRangeCompare comp = CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
     332  if (comp != CbcRangeOverlap || (comp==CbcRangeOverlap && !replaceIfOverlap)) {
     333    return comp;
     334  }
     335  r0.setLb(thisBd[0]);
     336  r0.setUb(thisBd[1]);
     337  return comp;
     338}
     339
     340//##############################################################################
    287341
    288342/** Default Constructor
  • trunk/Cbc/src/CbcBranchCut.hpp

    r765 r912  
    139139  virtual double branch();
    140140
     141#if 0
     142  // No need to override. Default works fine.
     143  /** Reset every information so that the branching object appears to point to
     144      the previous child. This method does not need to modify anything in any
     145      solver. */
     146  virtual void previousBranch();
     147#endif
     148
    141149  using CbcBranchingObject::print ;
    142150  /** \brief Print something about branch - only if log level high
     
    147155  */
    148156  virtual bool boundBranch() const;
     157
     158  /** Return the type (an integer identifier) of \c this */
     159  virtual int type() const { return 200; }
     160
     161  /** Compare the original object of \c this with the original object of \c
     162      brObj. Assumes that there is an ordering of the original objects.
     163      This method should be invoked only if \c this and brObj are of the same
     164      type.
     165      Return negative/0/positive depending on whether \c this is
     166      smaller/same/larger than the argument.
     167  */
     168  virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
     169
     170  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     171      same type and must have the same original object, but they may have
     172      different feasible regions.
     173      Return the appropriate CbcRangeCompare value (first argument being the
     174      sub/superset if that's the case). In case of overlap (and if \c
     175      replaceIfOverlap is true) replace the current branching object with one
     176      whose feasible region is the overlap.
     177   */
     178  virtual CbcRangeCompare compareBranchingObject
     179  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
    149180
    150181protected:
  • trunk/Cbc/src/CbcBranchDynamic.hpp

    r848 r912  
    393393  inline void setObject(CbcSimpleIntegerDynamicPseudoCost * object)
    394394  { object_=object;}
     395
     396  /** Return the type (an integer identifier) of \c this */
     397  virtual int type() const { return 400; }
     398
     399  // LL: compareOriginalObject and compareBranchingObject are inherited from
     400  // CbcIntegerBranchingObject thus need not be declared/defined here. After
     401  // all, this kind of branching object is simply using pseudocosts to make
     402  // decisions, but once the decisions are made they are the same kind as in
     403  // the underlying class.
     404
    395405protected:
    396406  /// Change in guessed objective value for next branch
  • trunk/Cbc/src/CbcBranchLotsize.cpp

    r904 r912  
    779779  }
    780780}
     781
     782/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     783    same type and must have the same original object, but they may have
     784    different feasible regions.
     785    Return the appropriate CbcRangeCompare value (first argument being the
     786    sub/superset if that's the case). In case of overlap (and if \c
     787    replaceIfOverlap is true) replace the current branching object with one
     788    whose feasible region is the overlap.
     789*/
     790CbcRangeCompare
     791CbcLotsizeBranchingObject::compareBranchingObject
     792(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
     793{
     794  const CbcLotsizeBranchingObject* br =
     795    dynamic_cast<const CbcLotsizeBranchingObject*>(brObj);
     796  assert(br);
     797  double* thisBd = way_ == -1 ? down_ : up_;
     798  const double* otherBd = br->way_ == -1 ? br->down_ : br->up_;
     799  return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
     800}
  • trunk/Cbc/src/CbcBranchLotsize.hpp

    r765 r912  
    197197  virtual double branch();
    198198
     199#if 0
     200  // No need to override. Default works fine.
     201  /** Reset every information so that the branching object appears to point to
     202      the previous child. This method does not need to modify anything in any
     203      solver. */
     204  virtual void previousBranch();
     205#endif
     206
    199207  using CbcBranchingObject::print ;
    200208  /** \brief Print something about branch - only if log level high
    201209  */
    202210  virtual void print();
     211
     212  /** Return the type (an integer identifier) of \c this */
     213  virtual int type() const { return 300; }
     214
     215  // LL: compareOriginalObject can be inherited from the CbcBranchingObject
     216  // since variable_ uniquely defines the lot sizing object.
     217
     218  /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
     219      same type and must have the same original object, but they may have
     220      different feasible regions.
     221      Return the appropriate CbcRangeCompare value (first argument being the
     222      sub/superset if that's the case). In case of overlap (and if \c
     223      replaceIfOverlap is true) replace the current branching object with one
     224      whose feasible region is the overlap.
     225   */
     226  virtual CbcRangeCompare compareBranchingObject
     227  (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
    203228
    204229protected:
  • trunk/Cbc/src/CbcHeuristic.cpp

    r904 r912  
    1313#include <cfloat>
    1414
     15//#define PRINT_DEBUG
    1516#ifdef COIN_HAS_CLP
    1617#include "OsiClpSolverInterface.hpp"
     
    2425#include "OsiAuxInfo.hpp"
    2526#include "OsiPresolve.hpp"
     27#include "CbcBranchActual.hpp"
     28//==============================================================================
     29
     30CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs)
     31{
     32  numObjects_ = rhs.numObjects_;
     33  brObj_ = new CbcBranchingObject*[numObjects_];
     34  for (int i = 0; i < numObjects_; ++i) {
     35    brObj_[i] = rhs.brObj_[i]->clone();
     36  }
     37}
     38
     39void
     40CbcHeuristicNodeList::gutsOfDelete()
     41{
     42  for (int i = nodes_.size() - 1; i >= 0; --i) {
     43    delete nodes_[i];
     44  }
     45}
     46
     47void
     48CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs)
     49{
     50  append(rhs);
     51}
     52
     53CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs)
     54{
     55  gutsOfCopy(rhs);
     56}
     57
     58CbcHeuristicNodeList& CbcHeuristicNodeList::operator=
     59(const CbcHeuristicNodeList& rhs)
     60{
     61  if (this != &rhs) {
     62    gutsOfDelete();
     63    gutsOfCopy(rhs);
     64  }
     65  return *this;
     66}
     67
     68CbcHeuristicNodeList::~CbcHeuristicNodeList()
     69{
     70  gutsOfDelete();
     71}
     72
     73void
     74CbcHeuristicNodeList::append(CbcHeuristicNode*& node)
     75{
     76  nodes_.push_back(node);
     77  node = NULL;
     78}
     79
     80void
     81CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes)
     82{
     83  nodes_.reserve(nodes_.size() + nodes.size());
     84  for (int i = 0; i < nodes.size(); ++i) {
     85    CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i));
     86    append(node);
     87  }
     88}
     89
     90//==============================================================================
    2691
    2792// Default Constructor
    28 CbcHeuristic::CbcHeuristic()
    29   :model_(NULL),
    30    when_(2),
    31    numberNodes_(200),
    32    feasibilityPumpOptions_(-1),
    33    fractionSmall_(1.0),
    34    heuristicName_("Unknown")
     93CbcHeuristic::CbcHeuristic() :
     94  model_(NULL),
     95  when_(2),
     96  numberNodes_(200),
     97  feasibilityPumpOptions_(-1),
     98  fractionSmall_(1.0),
     99  heuristicName_("Unknown"),
     100  howOften_(1),
     101  decayFactor_(0.0),
     102  shallowDepth_(1),
     103  howOftenShallow_(1),
     104  numInvocationsInShallow_(0),
     105  numInvocationsInDeep_(0),
     106  lastRunDeep_(0),
     107  numRuns_(0),
     108  minDistanceToRun_(1),
     109  runNodes_(),
     110  numCouldRun_(0)
    35111{
    36112  // As CbcHeuristic virtual need to modify .cpp if above change
     
    38114
    39115// Constructor from model
    40 CbcHeuristic::CbcHeuristic(CbcModel & model)
    41 :
     116CbcHeuristic::CbcHeuristic(CbcModel & model) :
    42117  model_(&model),
    43118  when_(2),
     
    45120  feasibilityPumpOptions_(-1),
    46121  fractionSmall_(1.0),
    47   heuristicName_("Unknown")
    48 {
    49   // As CbcHeuristic virtual need to modify .cpp if above change
     122  heuristicName_("Unknown"),
     123  howOften_(1),
     124  decayFactor_(0.0),
     125  shallowDepth_(1),
     126  howOftenShallow_(1),
     127  numInvocationsInShallow_(0),
     128  numInvocationsInDeep_(0),
     129  lastRunDeep_(0),
     130  numRuns_(0),
     131  minDistanceToRun_(1),
     132  runNodes_(),
     133  numCouldRun_(0)
     134{}
     135
     136void
     137CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs)
     138{
     139  model_ = rhs.model_;
     140  when_ = rhs.when_;
     141  numberNodes_ = rhs.numberNodes_;
     142  feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
     143  fractionSmall_ = rhs.fractionSmall_;
     144  randomNumberGenerator_ = rhs.randomNumberGenerator_;
     145  heuristicName_ = rhs.heuristicName_;
     146  howOften_ = rhs.howOften_;
     147  decayFactor_ = rhs.howOften_;
     148  shallowDepth_= rhs.shallowDepth_;
     149  howOftenShallow_= rhs.howOftenShallow_;
     150  numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
     151  numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
     152  lastRunDeep_ = rhs.lastRunDeep_;
     153  numRuns_ = rhs.numRuns_;
     154  numCouldRun_ = rhs.numCouldRun_;
     155  minDistanceToRun_ = rhs.minDistanceToRun_;
     156  runNodes_ = rhs.runNodes_;
    50157}
    51158// Copy constructor
    52159CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
    53 :
    54   model_(rhs.model_),
    55   when_(rhs.when_),
    56   numberNodes_(rhs.numberNodes_),
    57   feasibilityPumpOptions_(rhs.feasibilityPumpOptions_),
    58   fractionSmall_(rhs.fractionSmall_),
    59   randomNumberGenerator_(rhs.randomNumberGenerator_),
    60   heuristicName_(rhs.heuristicName_)
    61 {
    62 }
     160{
     161  gutsOfCopy(rhs);
     162}
     163
    63164// Assignment operator
    64165CbcHeuristic &
     
    66167{
    67168  if (this!=&rhs) {
    68     model_ = rhs.model_;
    69     when_ = rhs.when_;
    70     numberNodes_ = rhs.numberNodes_;
    71     feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
    72     fractionSmall_ = rhs.fractionSmall_;
    73     randomNumberGenerator_ = rhs.randomNumberGenerator_;
    74     heuristicName_ = rhs.heuristicName_ ;
     169    gutsOfDelete();
     170    gutsOfCopy(rhs);
    75171  }
    76172  return *this;
     173}
     174
     175void CbcHeurDebugNodes(CbcModel* model_)
     176{
     177  CbcNode* node = model_->currentNode();
     178  CbcNodeInfo* nodeInfo = node->nodeInfo();
     179  std::cout << "===============================================================\n";
     180  while (nodeInfo) {
     181    const CbcNode* node = nodeInfo->owner();
     182    printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
     183    {
     184      const CbcIntegerBranchingObject* brPrint =
     185        dynamic_cast<const CbcIntegerBranchingObject*>(nodeInfo->parentBranch());
     186      if (!brPrint) {
     187        printf("    parentBranch: NULL\n");
     188      } else {
     189        const double* downBounds = brPrint->downBounds();
     190        const double* upBounds = brPrint->upBounds();
     191        int variable = brPrint->variable();
     192        int way = brPrint->way();
     193        printf("   parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     194               variable, (int)downBounds[0], (int)downBounds[1],
     195               (int)upBounds[0], (int)upBounds[1], way);
     196      }
     197    }
     198    if (! node) {
     199      printf("    owner: NULL\n");
     200    } else {
     201      printf("    owner: node %i depth %i onTree %i active %i",
     202             node->nodeNumber(), node->depth(), node->onTree(), node->active());
     203      const OsiBranchingObject* osibr =
     204        nodeInfo->owner()->branchingObject();
     205      const CbcBranchingObject* cbcbr =
     206        dynamic_cast<const CbcBranchingObject*>(osibr);
     207      const CbcIntegerBranchingObject* brPrint =
     208        dynamic_cast<const CbcIntegerBranchingObject*>(cbcbr);
     209      if (!brPrint) {
     210        printf("        ownerBranch: NULL\n");
     211      } else {
     212        const double* downBounds = brPrint->downBounds();
     213        const double* upBounds = brPrint->upBounds();
     214        int variable = brPrint->variable();
     215        int way = brPrint->way();
     216        printf("        ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     217               variable, (int)downBounds[0], (int)downBounds[1],
     218               (int)upBounds[0], (int)upBounds[1], way);
     219      }
     220    }
     221    nodeInfo = nodeInfo->parent();
     222  }
     223}
     224
     225void
     226CbcHeuristic::debugNodes()
     227{
     228  CbcHeurDebugNodes(model_);
     229}
     230
     231void
     232CbcHeuristic::printDistanceToNodes()
     233{
     234  const CbcNode* currentNode = model_->currentNode();
     235  if (currentNode != NULL) {
     236    CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
     237    for (int i = runNodes_.size() - 1; i >= 0; --i) {
     238      nodeDesc->distance(runNodes_.node(i));
     239    }
     240    runNodes_.append(nodeDesc);
     241  }
     242}
     243
     244bool
     245CbcHeuristic::shouldHeurRun()
     246{
     247
     248#if 0
     249  const CbcNode* currentNode = model_->currentNode();
     250  if (currentNode == NULL) {
     251    return false;
     252  }
     253
     254  debugNodes();
     255//   return false;
     256
     257  const int depth = currentNode->depth();
     258#else
     259  int depth = 0;
     260  const CbcNode* currentNode = model_->currentNode();
     261  if (currentNode != NULL) {
     262    depth = currentNode->depth();
     263#ifdef PRINT_DEBUG
     264    debugNodes();
     265#endif
     266  }
     267#endif
     268
     269  const int nodeCount = model_->getNodeCount();  // FIXME: check that this is
     270                                                 // correct in parallel
     271
     272  if (nodeCount == 0 || depth <= shallowDepth_) {
     273    // what to do when we are in the shallow part of the tree
     274    if (model_->getCurrentPassNumber() == 1) {
     275      // first time in the node...
     276      numInvocationsInShallow_ = 0;
     277    }
     278    ++numInvocationsInShallow_;
     279    // Very large howOftenShallow_ will give the original test:
     280    // (model_->getCurrentPassNumber() != 1)
     281    //    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
     282    if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
     283      return false;
     284    }
     285    // LL: should we save these nodes in the list of nodes where the heur was
     286    // LL: run?
     287#if 1
     288    if (currentNode != NULL) {
     289      // Get where we are and create the appropriate CbcHeuristicNode object
     290      CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
     291      runNodes_.append(nodeDesc);
     292    }
     293#endif
     294  } else {
     295    // deeper in the tree
     296    if (model_->getCurrentPassNumber() == 1) {
     297      // first time in the node...
     298      ++numInvocationsInDeep_;
     299    }
     300    if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
     301      return false;
     302    }
     303    if (model_->getCurrentPassNumber() != 1) {
     304      // Run the heuristic only when first entering the node.
     305      // LL: I don't think this is right. It should run just before strong
     306      // LL: branching, I believe.
     307      return false;
     308    }
     309    // Get where we are and create the appropriate CbcHeuristicNode object
     310    CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
     311    //#ifdef PRINT_DEBUG
     312#if 1
     313    double minDistance = nodeDesc->minDistance(runNodes_);
     314    double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
     315    std::cout<<"minDistance = "<<minDistance
     316             <<", minDistanceToRun = "<<minDistanceToRun<<std::endl;
     317    if (minDistance < minDistanceToRun) {
     318#else
     319    if (nodeDesc->minDistance(runNodes_) < minDistanceToRun_) {
     320#endif
     321      delete nodeDesc;
     322      return false;
     323    }
     324    runNodes_.append(nodeDesc);
     325    lastRunDeep_ = numInvocationsInDeep_;
     326    //    ++lastRunDeep_;
     327  }
     328  ++numRuns_;
     329  return true;
     330}
     331
     332bool
     333CbcHeuristic::shouldHeurRun_randomChoice()
     334{
     335  int depth = 0;
     336  const CbcNode* currentNode = model_->currentNode();
     337  if (currentNode != NULL) {
     338    depth = currentNode->depth();
     339  }
     340
     341  if(depth != 0) {
     342    const double numerator = depth * depth;
     343    const double denominator = exp(depth * log((double)2));
     344    double probability = numerator / denominator;
     345    double randomNumber = randomNumberGenerator_.randomDouble();
     346    if (randomNumber>probability)
     347      return false;
     348   
     349    if (model_->getCurrentPassNumber() > 1)
     350      return false;
     351  }
     352  ++numRuns_;
     353  return true;
    77354}
    78355
     
    82359{
    83360  model_=model;
    84 }
     361  }
    85362// Set seed
    86363void
     
    460737  return returnCode;
    461738}
     739
     740//##############################################################################
     741
     742inline int compare3BranchingObjects(const CbcBranchingObject* br0,
     743                                    const CbcBranchingObject* br1)
     744{
     745  const int t0 = br0->type();
     746  const int t1 = br1->type();
     747  if (t0 < t1) {
     748    return -1;
     749  }
     750  if (t0 > t1) {
     751    return 1;
     752  }
     753  return br0->compareOriginalObject(br1);
     754}
     755
     756//==============================================================================
     757
     758inline bool compareBranchingObjects(const CbcBranchingObject* br0,
     759                                    const CbcBranchingObject* br1)
     760{
     761  return compare3BranchingObjects(br0, br1) < 0;
     762}
     763
     764//==============================================================================
     765
     766void
     767CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
     768{
     769  //  CbcHeurDebugNodes(&model);
     770  CbcNode* node = model.currentNode();
     771  brObj_ = new CbcBranchingObject*[node->depth()];
     772  CbcNodeInfo* nodeInfo = node->nodeInfo();
     773  int cnt = 0;
     774  while (nodeInfo->parentBranch() != NULL) {
     775    brObj_[cnt++] = nodeInfo->parentBranch()->clone();
     776    nodeInfo = nodeInfo->parent();
     777  }
     778  std::sort(brObj_, brObj_+cnt, compareBranchingObjects);
     779  if (cnt <= 1) {
     780    numObjects_ = cnt;
     781  } else {
     782    numObjects_ = 0;
     783    CbcBranchingObject* br;
     784    for (int i = 1; i < cnt; ++i) {
     785      if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
     786        int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], &br);
     787        switch (comp) {
     788        case CbcRangeSame: // the same range
     789        case CbcRangeDisjoint: // disjoint decisions
     790          // should not happen! we are on a chain!
     791          abort();
     792        case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
     793          delete brObj_[i];
     794          break;
     795        case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
     796          delete brObj_[numObjects_];
     797          brObj_[numObjects_] = brObj_[i];
     798          break;
     799        case CbcRangeOverlap: // overlap
     800          delete brObj_[i];
     801          delete brObj_[numObjects_];
     802          brObj_[numObjects_] = br;
     803          break;
     804      }
     805        continue;
     806      } else {
     807        brObj_[++numObjects_] = brObj_[i];
     808      }
     809    }
     810    ++numObjects_;
     811  }
     812}
     813
     814//==============================================================================
     815
     816CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
     817{
     818  gutsOfConstructor(model);
     819}
     820
     821//==============================================================================
     822
     823double
     824CbcHeuristicNode::distance(const CbcHeuristicNode* node) const
     825{
     826
     827  const double disjointWeight = 1;
     828  const double overlapWeight = 0.4;
     829  const double subsetWeight = 0.2;
     830  int countDisjointWeight = 0;
     831  int countOverlapWeight = 0;
     832  int countSubsetWeight = 0;
     833  int i = 0;
     834  int j = 0;
     835  double dist = 0.0;
     836#ifdef PRINT_DEBUG
     837  printf(" numObjects_ = %i, node->numObjects_ = %i\n",
     838         numObjects_, node->numObjects_);
     839#endif
     840  while( i < numObjects_ && j < node->numObjects_) {
     841    CbcBranchingObject* br0 = brObj_[i];
     842    const CbcBranchingObject* br1 = node->brObj_[j];
     843#ifdef PRINT_DEBUG
     844    const CbcIntegerBranchingObject* brPrint0 =
     845      dynamic_cast<const CbcIntegerBranchingObject*>(br0);
     846    const double* downBounds = brPrint0->downBounds();
     847    const double* upBounds = brPrint0->upBounds();
     848    int variable = brPrint0->variable();
     849    int way = brPrint0->way();
     850    printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     851           variable, (int)downBounds[0], (int)downBounds[1],
     852           (int)upBounds[0], (int)upBounds[1], way);
     853    const CbcIntegerBranchingObject* brPrint1 =
     854      dynamic_cast<const CbcIntegerBranchingObject*>(br1);
     855    downBounds = brPrint1->downBounds();
     856    upBounds = brPrint1->upBounds();
     857    variable = brPrint1->variable();
     858    way = brPrint1->way();
     859    printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     860           variable, (int)downBounds[0], (int)downBounds[1],
     861           (int)upBounds[0], (int)upBounds[1], way);
     862#endif
     863    const int brComp = compare3BranchingObjects(br0, br1);
     864    if (brComp < 0) {
     865      dist += subsetWeight;
     866      countSubsetWeight++;
     867      ++i;
     868    }
     869    else if (brComp > 0) {
     870      dist += subsetWeight;
     871      countSubsetWeight++;
     872      ++j;
     873    }
     874    else {
     875      const int comp = br0->compareBranchingObject(br1, false);
     876      switch (comp) {
     877      case CbcRangeSame:
     878        // do nothing
     879        break;
     880      case CbcRangeDisjoint: // disjoint decisions
     881        dist += disjointWeight;
     882        countDisjointWeight++;
     883        break;
     884      case CbcRangeSubset: // subset one way or another
     885      case CbcRangeSuperset:
     886        dist += subsetWeight;
     887        countSubsetWeight++;
     888        break;
     889      case CbcRangeOverlap: // overlap
     890        dist += overlapWeight;
     891        countOverlapWeight++;
     892        break;
     893      }
     894      ++i;
     895      ++j;
     896    }
     897  }
     898  dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
     899  countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
     900  printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
     901         countOverlapWeight, countDisjointWeight);
     902  return dist;
     903}
     904
     905//==============================================================================
     906
     907CbcHeuristicNode::~CbcHeuristicNode()
     908{
     909  for (int i = 0; i < numObjects_; ++i) {
     910    delete brObj_[i];
     911  }
     912  delete [] brObj_;
     913}
     914
     915//==============================================================================
     916
     917double
     918CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList)
     919{
     920  double minDist = COIN_DBL_MAX;
     921  for (int i = nodeList.size() - 1; i >= 0; --i) {
     922    minDist = CoinMin(minDist, distance(nodeList.node(i)));
     923  }
     924  return minDist;
     925}
     926
     927//==============================================================================
     928
     929double
     930CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList)
     931{
     932  if (nodeList.size() == 0) {
     933    return COIN_DBL_MAX;
     934  }
     935  double sumDist = 0;
     936  for (int i = nodeList.size() - 1; i >= 0; --i) {
     937    sumDist += distance(nodeList.node(i));
     938  }
     939  return sumDist/nodeList.size();
     940}
     941
     942//##############################################################################
    462943
    463944// Default Constructor
  • trunk/Cbc/src/CbcHeuristic.hpp

    r854 r912  
    99#include "OsiCuts.hpp"
    1010#include "CoinHelperFunctions.hpp"
     11#include "OsiBranchingObject.hpp"
    1112
    1213class OsiSolverInterface;
    1314
    1415class CbcModel;
     16
     17//#############################################################################
     18
     19class CbcHeuristicNodeList;
     20class CbcBranchingObject;
     21
     22/** A class describing the branching decisions that were made to get
     23    to the node where a heuristics was invoked from */
     24
     25class CbcHeuristicNode {
     26private:
     27  void gutsOfConstructor(CbcModel& model);
     28  CbcHeuristicNode();
     29  CbcHeuristicNode& operator=(const CbcHeuristicNode&);
     30private:
     31  /// The number of branching decisions made
     32  int numObjects_;
     33  /** The indices of the branching objects. Note: an index may be
     34      listed multiple times. E.g., a general integer variable that has
     35      been branched on multiple times. */
     36  CbcBranchingObject** brObj_;
     37public:
     38  CbcHeuristicNode(CbcModel& model);
     39
     40  CbcHeuristicNode(const CbcHeuristicNode& rhs);
     41  ~CbcHeuristicNode();
     42  double distance(const CbcHeuristicNode* node) const;
     43  double minDistance(const CbcHeuristicNodeList& nodeList);
     44  double avgDistance(const CbcHeuristicNodeList& nodeList);
     45};
     46
     47class CbcHeuristicNodeList {
     48private:
     49  void gutsOfDelete();
     50  void gutsOfCopy(const CbcHeuristicNodeList& rhs);
     51private:
     52  std::vector<CbcHeuristicNode*> nodes_;
     53public:
     54  CbcHeuristicNodeList() {}
     55  CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs);
     56  CbcHeuristicNodeList& operator=(const CbcHeuristicNodeList& rhs);
     57  ~CbcHeuristicNodeList();
     58 
     59  void append(CbcHeuristicNode*& node);
     60  void append(const CbcHeuristicNodeList& nodes);
     61  inline const CbcHeuristicNode* node(int i) const { return nodes_[i]; }
     62  inline int size() const { return nodes_.size(); }
     63};
    1564
    1665//#############################################################################
     
    1867
    1968class CbcHeuristic {
     69private:
     70  void gutsOfDelete() {}
     71  void gutsOfCopy(const CbcHeuristic & rhs);
     72
    2073public:
    2174  // Default Constructor
     
    124177  void setSeed(int value);
    125178
     179  /** Check whether the heuristic should run */
     180  bool shouldHeurRun();
     181  bool shouldHeurRun_randomChoice();
     182  void debugNodes();
     183  void printDistanceToNodes();
     184
    126185protected:
    127186
     
    140199  /// Name for printing
    141200  std::string heuristicName_;
    142  
     201
     202  /// How often to do (code can change)
     203  int howOften_;
     204  /// How much to increase how often
     205  double decayFactor_;
     206  /** Upto this depth we call the tree shallow and the heuristic can be called
     207      multiple times. That is, the test whether the current node is far from
     208      the others where the jeuristic was invoked will not be done, only the
     209      frequency will be tested. After that depth the heuristic will can be
     210      invoked only once per node, right before branching. That's when it'll be
     211      tested whether the heur should run at all. */
     212  int shallowDepth_;
     213  /** How often to invoke the heuristics in the shallow part of the tree */
     214  int howOftenShallow_;
     215  /** How many invocations happened within the same node when in a shallow
     216      part of the tree. */
     217  int numInvocationsInShallow_;
     218  /** How many invocations happened when in the deep part of the tree. For
     219      every node we count only one invocation. */
     220  int numInvocationsInDeep_;
     221  /** After how many deep invocations was the heuristic run last time */
     222  int lastRunDeep_;
     223  /// how many times the heuristic has actually run
     224  int numRuns_;
     225  /** How "far" should this node be from every other where the heuristic was
     226      run in order to allow the heuristic to run in this node, too. Currently
     227      this is tested, but we may switch to avgDistanceToRun_ in the future. */
     228  int minDistanceToRun_;
     229
     230  /// The description of the nodes where this heuristic has been applied
     231  CbcHeuristicNodeList runNodes_;
     232
     233  /// How many times the heuristic could run
     234  int numCouldRun_;
     235
     236
     237#if 0
     238  /// Lower bounds of last node where the heuristic found a solution
     239  double * lowerBoundLastNode_;
     240  /// Upper bounds of last node where the heuristic found a solution
     241  double * upperBoundLastNode_;
     242#endif
    143243};
    144244/** Rounding class
  • trunk/Cbc/src/CbcHeuristicDiveCoefficient.cpp

    r873 r912  
    66#endif
    77
     8//#define PRINT_DEBUG
     9
    810#include "CbcHeuristicDiveCoefficient.hpp"
    911#include "CbcStrategy.hpp"
    10 #include  "CoinTime.hpp"
    1112
    1213// Default Constructor
    1314CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient()
    14   :CbcHeuristic()
     15  :CbcHeuristicDive()
    1516{
    16   // matrix and row copy will automatically be empty
    17   downLocks_ =NULL;
    18   upLocks_ = NULL;
    19   percentageToFix_ = 0.2;
    20   maxIterations_ = 100;
    21   maxTime_ = 60;
    2217}
    2318
    2419// Constructor from model
    2520CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(CbcModel & model)
    26   :CbcHeuristic(model)
     21  :CbcHeuristicDive(model)
    2722{
    28   downLocks_ =NULL;
    29   upLocks_ = NULL;
    30   // Get a copy of original matrix
    31   assert(model.solver());
    32   // model may have empty matrix - wait until setModel
    33   const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
    34   if (matrix) {
    35     matrix_ = *matrix;
    36   }
    37   percentageToFix_ = 0.2;
    38   maxIterations_ = 100;
    39   maxTime_ = 60;
    4023}
    4124
     
    4326CbcHeuristicDiveCoefficient::~CbcHeuristicDiveCoefficient ()
    4427{
    45   delete [] downLocks_;
    46   delete [] upLocks_;
    4728}
    4829
     
    6243  fprintf(fp,"3  CbcHeuristicDiveCoefficient heuristicDiveCoefficient(*cbcModel);\n");
    6344  CbcHeuristic::generateCpp(fp,"heuristicDiveCoefficient");
    64   if (percentageToFix_!=other.percentageToFix_)
    65     fprintf(fp,"3  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
    66   else
    67     fprintf(fp,"4  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
    68   if (maxIterations_!=other.maxIterations_)
    69     fprintf(fp,"3  heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);
    70   else
    71     fprintf(fp,"4  heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);
    72   if (maxTime_!=other.maxTime_)
    73     fprintf(fp,"3  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
    74   else
    75     fprintf(fp,"4  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
    7645  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveCoefficient);\n");
    7746}
     
    8049CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(const CbcHeuristicDiveCoefficient & rhs)
    8150:
    82   CbcHeuristic(rhs),
    83   matrix_(rhs.matrix_),
    84   percentageToFix_(rhs.percentageToFix_),
    85   maxIterations_(rhs.maxIterations_),
    86   maxTime_(rhs.maxTime_)
     51  CbcHeuristicDive(rhs)
    8752{
    88   int numberIntegers = model_->numberIntegers();
    89   if (rhs.downLocks_) {
    90     int numberIntegers = model_->numberIntegers();
    91     downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    92     upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    93   } else {
    94     downLocks_ = NULL;
    95     upLocks_ = NULL;
    96   }
    9753}
    9854
     
    10258{
    10359  if (this!=&rhs) {
    104     CbcHeuristic::operator=(rhs);
    105     matrix_ = rhs.matrix_;
    106     percentageToFix_ = rhs.percentageToFix_;
    107     maxIterations_ = rhs.maxIterations_;
    108     maxTime_ = rhs.maxTime_;
    109     delete [] downLocks_;
    110     delete [] upLocks_;
    111     if (rhs.downLocks_) {
    112       int numberIntegers = model_->numberIntegers();
    113       downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    114       upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    115     } else {
    116       downLocks_ = NULL;
    117       upLocks_ = NULL;
    118     }
     60    CbcHeuristicDive::operator=(rhs);
    11961  }
    12062  return *this;
    12163}
    12264
    123 // Resets stuff if model changes
    124 void
    125 CbcHeuristicDiveCoefficient::resetModel(CbcModel * model)
     65void
     66CbcHeuristicDiveCoefficient::selectVariableToBranch(OsiSolverInterface* solver,
     67                                                    const double* newSolution,
     68                                                    int& bestColumn,
     69                                                    int& bestRound)
    12670{
    127   model_=model;
    128   // Get a copy of original matrix
    129   assert(model_->solver());
    130   // model may have empty matrix - wait until setModel
    131   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    132   if (matrix) {
    133     matrix_ = *matrix;
    134     validate();
    135   }
    136 }
    137 
    138 // See if dive fractional will give better solution
    139 // Sets value of solution
    140 // Returns 1 if solution, 0 if not
    141 int
    142 CbcHeuristicDiveCoefficient::solution(double & solutionValue,
    143                                      double * betterSolution)
    144 {
    145 
    146   // See if to do
    147   if (!when()||(when()%10==1&&model_->phase()!=1)||
    148       (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
    149     return 0; // switched off
    150 
    151   double time1 = CoinCpuTime();
    152 
    153   OsiSolverInterface * solver = model_->solver()->clone();
    154   const double * lower = solver->getColLower();
    155   const double * upper = solver->getColUpper();
    156   const double * rowLower = solver->getRowLower();
    157   const double * rowUpper = solver->getRowUpper();
    158   const double * solution = solver->getColSolution();
    159   const double * objective = solver->getObjCoefficients();
    160   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    161   double primalTolerance;
    162   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    163 
    164   int numberRows = matrix_.getNumRows();
    165   assert (numberRows<=solver->getNumRows());
    16671  int numberIntegers = model_->numberIntegers();
    16772  const int * integerVariable = model_->integerVariable();
    168   double direction = solver->getObjSense(); // 1 for min, -1 for max
    169   double newSolutionValue = direction*solver->getObjValue();
    170   int returnCode = 0;
    171   // Column copy
    172   const double * element = matrix_.getElements();
    173   const int * row = matrix_.getIndices();
    174   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    175   const int * columnLength = matrix_.getVectorLengths();
     73  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    17674
    177   // Get solution array for heuristic solution
    178   int numberColumns = solver->getNumCols();
    179   double * newSolution = new double [numberColumns];
    180   memcpy(newSolution,solution,numberColumns*sizeof(double));
    181 
    182   // vectors to store the latest variables fixed at their bounds
    183   int* columnFixed = new int [numberIntegers];
    184   double* originalBound = new double [numberIntegers];
    185   bool * fixedAtLowerBound = new bool [numberIntegers];
    186 
    187   const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
    188 
    189   // count how many fractional variables
    190   int numberFractionalVariables = 0;
     75  bestColumn = -1;
     76  bestRound = -1; // -1 rounds down, +1 rounds up
     77  double bestFraction = DBL_MAX;
     78  int bestLocks = COIN_INT_MAX;
    19179  for (int i=0; i<numberIntegers; i++) {
    19280    int iColumn = integerVariable[i];
    19381    double value=newSolution[iColumn];
     82    double fraction=value-floor(value);
     83    int round = 0;
    19484    if (fabs(floor(value+0.5)-value)>integerTolerance) {
    195       numberFractionalVariables++;
    196     }
    197   }
    198 
    199   int iteration = 0;
    200   while(numberFractionalVariables) {
    201     iteration++;
    202 
    203     // select a fractional variable to bound
    204     double bestFraction = DBL_MAX;
    205     int bestColumn = -1;
    206     int bestRound = -1; // -1 rounds down, +1 rounds up
    207     int bestLocks = COIN_INT_MAX;
    208     int numberAtBoundFixed = 0;
    209     bool canRoundSolution = true;
    210     for (int i=0; i<numberIntegers; i++) {
    211       int iColumn = integerVariable[i];
    212       double value=newSolution[iColumn];
    213       double fraction=value-floor(value);
    214       int round = 0;
    215       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    216         int nDownLocks = downLocks_[i];
    217         int nUpLocks = upLocks_[i];
    218         if(nDownLocks>0 && nUpLocks>0) {
    219           // the variable cannot be rounded
    220           canRoundSolution = false;
    221           int nLocks = nDownLocks;
    222           if(nDownLocks<nUpLocks)
    223             round = -1;
    224           else if(nDownLocks>nUpLocks) {
    225             round = 1;
    226             fraction = 1.0 - fraction;
    227             nLocks = nUpLocks;
    228           }
    229           else if(fraction < 0.5)
    230             round = -1;
    231           else {
    232             round = 1;
    233             fraction = 1.0 - fraction;
    234             nLocks = nUpLocks;
    235           }
    236 
    237           // if variable is not binary, penalize it
    238           if(!solver->isBinary(iColumn))
    239             fraction *= 1000.0;
    240 
    241           if(nLocks < bestLocks || (nLocks == bestLocks &&
    242                                     fraction < bestFraction)) {
    243             bestColumn = iColumn;
    244             bestLocks = nLocks;
    245             bestFraction = fraction;
    246             bestRound = round;
    247           }
     85      int nDownLocks = downLocks_[i];
     86      int nUpLocks = upLocks_[i];
     87      if(nDownLocks>0 && nUpLocks>0) {
     88        // the variable cannot be rounded
     89        int nLocks = nDownLocks;
     90        if(nDownLocks<nUpLocks)
     91          round = -1;
     92        else if(nDownLocks>nUpLocks) {
     93          round = 1;
     94          fraction = 1.0 - fraction;
     95          nLocks = nUpLocks;
    24896        }
    249       }
    250       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    251         // fix the variable at one of its bounds
    252         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    253             lower[iColumn] != upper[iColumn]) {
    254           columnFixed[numberAtBoundFixed] = iColumn;
    255           originalBound[numberAtBoundFixed] = upper[iColumn];
    256           fixedAtLowerBound[numberAtBoundFixed] = true;
    257           solver->setColUpper(iColumn, lower[iColumn]);
    258           numberAtBoundFixed++;
     97        else if(fraction < 0.5)
     98          round = -1;
     99        else {
     100          round = 1;
     101          fraction = 1.0 - fraction;
     102          nLocks = nUpLocks;
    259103        }
    260         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    261             lower[iColumn] != upper[iColumn]) {
    262           columnFixed[numberAtBoundFixed] = iColumn;
    263           originalBound[numberAtBoundFixed] = lower[iColumn];
    264           fixedAtLowerBound[numberAtBoundFixed] = false;
    265           solver->setColLower(iColumn, upper[iColumn]);
    266           numberAtBoundFixed++;
     104       
     105        // if variable is not binary, penalize it
     106        if(!solver->isBinary(iColumn))
     107          fraction *= 1000.0;
     108       
     109        if(nLocks < bestLocks || (nLocks == bestLocks &&
     110                                  fraction < bestFraction)) {
     111          bestColumn = iColumn;
     112          bestLocks = nLocks;
     113          bestFraction = fraction;
     114          bestRound = round;
    267115        }
    268116      }
    269117    }
    270 
    271     if(canRoundSolution) {
    272       // Round all the fractional variables
    273       for (int i=0; i<numberIntegers; i++) {
    274         int iColumn = integerVariable[i];
    275         double value=newSolution[iColumn];
    276         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    277           if(downLocks_[i]==0 || upLocks_[i]==0) {
    278             if(downLocks_[i]==0 && upLocks_[i]==0) {
    279               if(direction * objective[iColumn] >= 0.0)
    280                 newSolution[iColumn] = floor(value);
    281               else
    282                 newSolution[iColumn] = ceil(value);
    283             }
    284             else if(downLocks_[i]==0)
    285               newSolution[iColumn] = floor(value);
    286             else
    287               newSolution[iColumn] = ceil(value);
    288           }
    289           else
    290             break;
    291         }
    292       }
    293       break;
    294     }
    295 
    296 
    297     double originalBoundBestColumn;
    298     if(bestColumn >= 0) {
    299       if(bestRound < 0) {
    300         originalBoundBestColumn = upper[bestColumn];
    301         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    302       }
    303       else {
    304         originalBoundBestColumn = lower[bestColumn];
    305         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    306       }
    307     } else {
    308       break;
    309     }
    310     int originalBestRound = bestRound;
    311     while (1) {
    312 
    313       solver->resolve();
    314 
    315       if(!solver->isProvenOptimal()) {
    316         if(numberAtBoundFixed > 0) {
    317           // Remove the bound fix for variables that were at bounds
    318           for(int i=0; i<numberAtBoundFixed; i++) {
    319             int iColFixed = columnFixed[i];
    320             if(fixedAtLowerBound[i])
    321               solver->setColUpper(iColFixed, originalBound[i]);
    322             else
    323               solver->setColLower(iColFixed, originalBound[i]);
    324           }
    325           numberAtBoundFixed = 0;
    326         }
    327         else if(bestRound == originalBestRound) {
    328           bestRound *= (-1);
    329           if(bestRound < 0) {
    330             solver->setColLower(bestColumn, originalBoundBestColumn);
    331             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    332           }
    333           else {
    334             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    335             solver->setColUpper(bestColumn, originalBoundBestColumn);
    336           }
    337         }
    338         else
    339           break;
    340       }
    341       else
    342         break;
    343     }
    344 
    345     if(!solver->isProvenOptimal())
    346       break;
    347 
    348     if(iteration > maxIterations_) {
    349       break;
    350     }
    351 
    352     if(CoinCpuTime()-time1 > maxTime_) {
    353       break;
    354     }
    355 
    356     memcpy(newSolution,solution,numberColumns*sizeof(double));
    357     numberFractionalVariables = 0;
    358     for (int i=0; i<numberIntegers; i++) {
    359       int iColumn = integerVariable[i];
    360       double value=newSolution[iColumn];
    361       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    362         numberFractionalVariables++;
    363       }
    364     }
    365 
    366   }
    367 
    368 
    369   double * rowActivity = new double[numberRows];
    370   memset(rowActivity,0,numberRows*sizeof(double));
    371 
    372   // re-compute new solution value
    373   double objOffset=0.0;
    374   solver->getDblParam(OsiObjOffset,objOffset);
    375   newSolutionValue = -objOffset;
    376   for (int i=0 ; i<numberColumns ; i++ )
    377     newSolutionValue += objective[i]*newSolution[i];
    378   newSolutionValue *= direction;
    379     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    380   if (newSolutionValue<solutionValue) {
    381     // paranoid check
    382     memset(rowActivity,0,numberRows*sizeof(double));
    383     for (int i=0;i<numberColumns;i++) {
    384       int j;
    385       double value = newSolution[i];
    386       if (value) {
    387         for (j=columnStart[i];
    388              j<columnStart[i]+columnLength[i];j++) {
    389           int iRow=row[j];
    390           rowActivity[iRow] += value*element[j];
    391         }
    392       }
    393     }
    394     // check was approximately feasible
    395     bool feasible=true;
    396     for (int i=0;i<numberRows;i++) {
    397       if(rowActivity[i]<rowLower[i]) {
    398         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    399           feasible = false;
    400       } else if(rowActivity[i]>rowUpper[i]) {
    401         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    402           feasible = false;
    403       }
    404     }
    405     for (int i=0; i<numberIntegers; i++) {
    406       int iColumn = integerVariable[i];
    407       double value=newSolution[iColumn];
    408       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    409         feasible = false;
    410         break;
    411       }
    412     }
    413     if (feasible) {
    414       // new solution
    415       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    416       solutionValue = newSolutionValue;
    417       //printf("** Solution of %g found by CbcHeuristicDiveCoefficient\n",newSolutionValue);
    418       returnCode=1;
    419     } else {
    420       // Can easily happen
    421       //printf("Debug CbcHeuristicDiveCoefficient giving bad solution\n");
    422     }
    423   }
    424 
    425   delete [] newSolution;
    426   delete [] columnFixed;
    427   delete [] originalBound;
    428   delete [] fixedAtLowerBound;
    429   delete [] rowActivity;
    430   delete solver;
    431   return returnCode;
    432 }
    433 
    434 // update model
    435 void CbcHeuristicDiveCoefficient::setModel(CbcModel * model)
    436 {
    437   model_ = model;
    438   // Get a copy of original matrix (and by row for rounding);
    439   assert(model_->solver());
    440   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    441   if (matrix) {
    442     matrix_ = *matrix;
    443     // make sure model okay for heuristic
    444     validate();
    445118  }
    446119}
    447 
    448 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    449 void
    450 CbcHeuristicDiveCoefficient::validate()
    451 {
    452   if (model_&&when()<10) {
    453     if (model_->numberIntegers()!=
    454         model_->numberObjects())
    455       setWhen(0);
    456   }
    457 
    458   int numberIntegers = model_->numberIntegers();
    459   const int * integerVariable = model_->integerVariable();
    460   delete [] downLocks_;
    461   delete [] upLocks_;
    462   downLocks_ = new unsigned short [numberIntegers];
    463   upLocks_ = new unsigned short [numberIntegers];
    464   // Column copy
    465   const double * element = matrix_.getElements();
    466   const int * row = matrix_.getIndices();
    467   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    468   const int * columnLength = matrix_.getVectorLengths();
    469   const double * rowLower = model_->solver()->getRowLower();
    470   const double * rowUpper = model_->solver()->getRowUpper();
    471   for (int i=0;i<numberIntegers;i++) {
    472     int iColumn = integerVariable[i];
    473     int down=0;
    474     int up=0;
    475     if (columnLength[iColumn]>65535) {
    476       setWhen(0);
    477       break; // unlikely to work
    478     }
    479     for (CoinBigIndex j=columnStart[iColumn];
    480          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    481       int iRow=row[j];
    482       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    483         up++;
    484         down++;
    485       } else if (element[j]>0.0) {
    486         if (rowUpper[iRow]<1.0e20)
    487           up++;
    488         else
    489           down++;
    490       } else {
    491         if (rowLower[iRow]>-1.0e20)
    492           up++;
    493         else
    494           down++;
    495       }
    496     }
    497     downLocks_[i] = (unsigned short) down;
    498     upLocks_[i] = (unsigned short) up;
    499   }
    500 }
  • trunk/Cbc/src/CbcHeuristicDiveCoefficient.hpp

    r868 r912  
    44#define CbcHeuristicDiveCoefficient_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveCoefficient class
    99 */
    1010
    11 class CbcHeuristicDiveCoefficient : public CbcHeuristic {
     11class CbcHeuristicDiveCoefficient : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
    37 
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does Coefficient Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     35  /// Selects the next variable to branch on
     36  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     37                                      const double* newSolution,
     38                                      int& bestColumn,
     39                                      int& bestRound);
    8640
    8741};
  • trunk/Cbc/src/CbcHeuristicDiveFractional.cpp

    r871 r912  
    88#include "CbcHeuristicDiveFractional.hpp"
    99#include "CbcStrategy.hpp"
    10 #include  "CoinTime.hpp"
    1110
    1211// Default Constructor
    1312CbcHeuristicDiveFractional::CbcHeuristicDiveFractional()
    14   :CbcHeuristic()
     13  :CbcHeuristicDive()
    1514{
    16   // matrix and row copy will automatically be empty
    17   downLocks_ =NULL;
    18   upLocks_ = NULL;
    19   percentageToFix_ = 0.2;
    20   maxIterations_ = 100;
    21   maxTime_ = 60;
    2215}
    2316
    2417// Constructor from model
    2518CbcHeuristicDiveFractional::CbcHeuristicDiveFractional(CbcModel & model)
    26   :CbcHeuristic(model)
     19  :CbcHeuristicDive(model)
    2720{
    28   downLocks_ =NULL;
    29   upLocks_ = NULL;
    30   // Get a copy of original matrix
    31   assert(model.solver());
    32   // model may have empty matrix - wait until setModel
    33   const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
    34   if (matrix) {
    35     matrix_ = *matrix;
    36   }
    37   percentageToFix_ = 0.2;
    38   maxIterations_ = 100;
    39   maxTime_ = 60;
    4021}
    4122
     
    4324CbcHeuristicDiveFractional::~CbcHeuristicDiveFractional ()
    4425{
    45   delete [] downLocks_;
    46   delete [] upLocks_;
    4726}
    4827
     
    6241  fprintf(fp,"3  CbcHeuristicDiveFractional heuristicDiveFractional(*cbcModel);\n");
    6342  CbcHeuristic::generateCpp(fp,"heuristicDiveFractional");
    64   if (percentageToFix_!=other.percentageToFix_)
    65     fprintf(fp,"3  heuristicDiveFractional.setPercentageToFix(%.f);\n",percentageToFix_);
    66   else
    67     fprintf(fp,"4  heuristicDiveFractional.setPercentageToFix(%.f);\n",percentageToFix_);
    68   if (maxIterations_!=other.maxIterations_)
    69     fprintf(fp,"3  heuristicDiveFractional.setMaxIterations(%d);\n",maxIterations_);
    70   else
    71     fprintf(fp,"4  heuristicDiveFractional.setMaxIterations(%d);\n",maxIterations_);
    72   if (maxTime_!=other.maxTime_)
    73     fprintf(fp,"3  heuristicDiveFractional.setMaxTime(%.2f);\n",maxTime_);
    74   else
    75     fprintf(fp,"4  heuristicDiveFractional.setMaxTime(%.2f);\n",maxTime_);
    7643  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveFractional);\n");
    7744}
     
    8047CbcHeuristicDiveFractional::CbcHeuristicDiveFractional(const CbcHeuristicDiveFractional & rhs)
    8148:
    82   CbcHeuristic(rhs),
    83   matrix_(rhs.matrix_),
    84   percentageToFix_(rhs.percentageToFix_),
    85   maxIterations_(rhs.maxIterations_),
    86   maxTime_(rhs.maxTime_)
     49  CbcHeuristicDive(rhs)
    8750{
    88   int numberIntegers = model_->numberIntegers();
    89   if (rhs.downLocks_) {
    90     int numberIntegers = model_->numberIntegers();
    91     downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    92     upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    93   } else {
    94     downLocks_ = NULL;
    95     upLocks_ = NULL;
    96   }
    9751}
    9852
     
    10256{
    10357  if (this!=&rhs) {
    104     CbcHeuristic::operator=(rhs);
    105     matrix_ = rhs.matrix_;
    106     percentageToFix_ = rhs.percentageToFix_;
    107     maxIterations_ = rhs.maxIterations_;
    108     maxTime_ = rhs.maxTime_;
    109     delete [] downLocks_;
    110     delete [] upLocks_;
    111     if (rhs.downLocks_) {
    112       int numberIntegers = model_->numberIntegers();
    113       downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    114       upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    115     } else {
    116       downLocks_ = NULL;
    117       upLocks_ = NULL;
    118     }
     58    CbcHeuristicDive::operator=(rhs);
    11959  }
    12060  return *this;
    12161}
    12262
    123 // Resets stuff if model changes
    124 void
    125 CbcHeuristicDiveFractional::resetModel(CbcModel * model)
     63void
     64CbcHeuristicDiveFractional::selectVariableToBranch(OsiSolverInterface* solver,
     65                                                   const double* newSolution,
     66                                                   int& bestColumn,
     67                                                   int& bestRound)
    12668{
    127   model_=model;
    128   // Get a copy of original matrix
    129   assert(model_->solver());
    130   // model may have empty matrix - wait until setModel
    131   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    132   if (matrix) {
    133     matrix_ = *matrix;
    134     validate();
    135   }
    136 }
    137 
    138 // See if dive fractional will give better solution
    139 // Sets value of solution
    140 // Returns 1 if solution, 0 if not
    141 int
    142 CbcHeuristicDiveFractional::solution(double & solutionValue,
    143                                      double * betterSolution)
    144 {
    145 
    146   // See if to do
    147   if (!when()||(when()%10==1&&model_->phase()!=1)||
    148       (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
    149     return 0; // switched off
    150 
    151   double time1 = CoinCpuTime();
    152 
    153   OsiSolverInterface * solver = model_->solver()->clone();
    154   const double * lower = solver->getColLower();
    155   const double * upper = solver->getColUpper();
    156   const double * rowLower = solver->getRowLower();
    157   const double * rowUpper = solver->getRowUpper();
    158   const double * solution = solver->getColSolution();
    159   const double * objective = solver->getObjCoefficients();
    160   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    161   double primalTolerance;
    162   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    163 
    164   int numberRows = matrix_.getNumRows();
    165   assert (numberRows<=solver->getNumRows());
    16669  int numberIntegers = model_->numberIntegers();
    16770  const int * integerVariable = model_->integerVariable();
    168   double direction = solver->getObjSense(); // 1 for min, -1 for max
    169   double newSolutionValue = direction*solver->getObjValue();
    170   int returnCode = 0;
    171   // Column copy
    172   const double * element = matrix_.getElements();
    173   const int * row = matrix_.getIndices();
    174   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    175   const int * columnLength = matrix_.getVectorLengths();
     71  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    17672
    177   // Get solution array for heuristic solution
    178   int numberColumns = solver->getNumCols();
    179   double * newSolution = new double [numberColumns];
    180   memcpy(newSolution,solution,numberColumns*sizeof(double));
    181 
    182   // vectors to store the latest variables fixed at their bounds
    183   int* columnFixed = new int [numberIntegers];
    184   double* originalBound = new double [numberIntegers];
    185   bool * fixedAtLowerBound = new bool [numberIntegers];
    186 
    187   const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
    188 
    189   // count how many fractional variables
    190   int numberFractionalVariables = 0;
     73  bestColumn = -1;
     74  bestRound = -1; // -1 rounds down, +1 rounds up
     75  double bestFraction = DBL_MAX;
    19176  for (int i=0; i<numberIntegers; i++) {
    19277    int iColumn = integerVariable[i];
    19378    double value=newSolution[iColumn];
     79    double fraction=value-floor(value);
     80    int round = 0;
    19481    if (fabs(floor(value+0.5)-value)>integerTolerance) {
    195       numberFractionalVariables++;
    196     }
    197   }
    198 
    199   int iteration = 0;
    200   while(numberFractionalVariables) {
    201     iteration++;
    202 
    203     // select a fractional variable to bound
    204     double bestFraction = DBL_MAX;
    205     int bestColumn = -1;
    206     int bestRound = -1; // -1 rounds down, +1 rounds up
    207     int numberAtBoundFixed = 0;
    208     bool canRoundSolution = true;
    209     for (int i=0; i<numberIntegers; i++) {
    210       int iColumn = integerVariable[i];
    211       double value=newSolution[iColumn];
    212       double fraction=value-floor(value);
    213       int round = 0;
    214       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    215         if(downLocks_[i]>0&&upLocks_[i]>0) {
     82      if(downLocks_[i]>0&&upLocks_[i]>0) {
    21683          // the variable cannot be rounded
    217           canRoundSolution = false;
    218           if(fraction < 0.5)
    219             round = -1;
    220           else {
    221             round = 1;
    222             fraction = 1.0 - fraction;
    223           }
    224 
    225           // if variable is not binary, penalize it
    226           if(!solver->isBinary(iColumn))
     84        if(fraction < 0.5)
     85          round = -1;
     86        else {
     87          round = 1;
     88          fraction = 1.0 - fraction;
     89        }
     90       
     91        // if variable is not binary, penalize it
     92        if(!solver->isBinary(iColumn))
    22793            fraction *= 1000.0;
    228 
    229           if(fraction < bestFraction) {
    230             bestColumn = iColumn;
    231             bestFraction = fraction;
    232             bestRound = round;
    233           }
    234         }
    235       }
    236       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    237         // fix the variable at one of its bounds
    238         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    239             lower[iColumn] != upper[iColumn]) {
    240           columnFixed[numberAtBoundFixed] = iColumn;
    241           originalBound[numberAtBoundFixed] = upper[iColumn];
    242           fixedAtLowerBound[numberAtBoundFixed] = true;
    243           solver->setColUpper(iColumn, lower[iColumn]);
    244           numberAtBoundFixed++;
    245         }
    246         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    247             lower[iColumn] != upper[iColumn]) {
    248           columnFixed[numberAtBoundFixed] = iColumn;
    249           originalBound[numberAtBoundFixed] = lower[iColumn];
    250           fixedAtLowerBound[numberAtBoundFixed] = false;
    251           solver->setColLower(iColumn, upper[iColumn]);
    252           numberAtBoundFixed++;
     94       
     95        if(fraction < bestFraction) {
     96          bestColumn = iColumn;
     97          bestFraction = fraction;
     98          bestRound = round;
    25399        }
    254100      }
    255101    }
    256 
    257     if(canRoundSolution) {
    258       // Round all the fractional variables
    259       for (int i=0; i<numberIntegers; i++) {
    260         int iColumn = integerVariable[i];
    261         double value=newSolution[iColumn];
    262         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    263           if(downLocks_[i]==0 || upLocks_[i]==0) {
    264             if(downLocks_[i]==0 && upLocks_[i]==0) {
    265               if(direction * objective[iColumn] >= 0.0)
    266                 newSolution[iColumn] = floor(value);
    267               else
    268                 newSolution[iColumn] = ceil(value);
    269             }
    270             else if(downLocks_[i]==0)
    271               newSolution[iColumn] = floor(value);
    272             else
    273               newSolution[iColumn] = ceil(value);
    274           }
    275           else
    276             break;
    277         }
    278       }
    279       break;
    280     }
    281 
    282 
    283     double originalBoundBestColumn;
    284     if(bestColumn >= 0) {
    285       if(bestRound < 0) {
    286         originalBoundBestColumn = upper[bestColumn];
    287         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    288       }
    289       else {
    290         originalBoundBestColumn = lower[bestColumn];
    291         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    292       }
    293     } else {
    294       break;
    295     }
    296     int originalBestRound = bestRound;
    297     while (1) {
    298 
    299       solver->resolve();
    300 
    301       if(!solver->isProvenOptimal()) {
    302         if(numberAtBoundFixed > 0) {
    303           // Remove the bound fix for variables that were at bounds
    304           for(int i=0; i<numberAtBoundFixed; i++) {
    305             int iColFixed = columnFixed[i];
    306             if(fixedAtLowerBound[i])
    307               solver->setColUpper(iColFixed, originalBound[i]);
    308             else
    309               solver->setColLower(iColFixed, originalBound[i]);
    310           }
    311           numberAtBoundFixed = 0;
    312         }
    313         else if(bestRound == originalBestRound) {
    314           bestRound *= (-1);
    315           if(bestRound < 0) {
    316             solver->setColLower(bestColumn, originalBoundBestColumn);
    317             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    318           }
    319           else {
    320             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    321             solver->setColUpper(bestColumn, originalBoundBestColumn);
    322           }
    323         }
    324         else
    325           break;
    326       }
    327       else
    328         break;
    329     }
    330 
    331     if(!solver->isProvenOptimal())
    332       break;
    333 
    334     if(iteration > maxIterations_) {
    335       break;
    336     }
    337 
    338     if(CoinCpuTime()-time1 > maxTime_) {
    339       break;
    340     }
    341 
    342     memcpy(newSolution,solution,numberColumns*sizeof(double));
    343     numberFractionalVariables = 0;
    344     for (int i=0; i<numberIntegers; i++) {
    345       int iColumn = integerVariable[i];
    346       double value=newSolution[iColumn];
    347       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    348         numberFractionalVariables++;
    349       }
    350     }
    351 
    352   }
    353 
    354 
    355   double * rowActivity = new double[numberRows];
    356   memset(rowActivity,0,numberRows*sizeof(double));
    357 
    358   // re-compute new solution value
    359   double objOffset=0.0;
    360   solver->getDblParam(OsiObjOffset,objOffset);
    361   newSolutionValue = -objOffset;
    362   for (int i=0 ; i<numberColumns ; i++ )
    363     newSolutionValue += objective[i]*newSolution[i];
    364   newSolutionValue *= direction;
    365     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    366   if (newSolutionValue<solutionValue) {
    367     // paranoid check
    368     memset(rowActivity,0,numberRows*sizeof(double));
    369     for (int i=0;i<numberColumns;i++) {
    370       int j;
    371       double value = newSolution[i];
    372       if (value) {
    373         for (j=columnStart[i];
    374              j<columnStart[i]+columnLength[i];j++) {
    375           int iRow=row[j];
    376           rowActivity[iRow] += value*element[j];
    377         }
    378       }
    379     }
    380     // check was approximately feasible
    381     bool feasible=true;
    382     for (int i=0;i<numberRows;i++) {
    383       if(rowActivity[i]<rowLower[i]) {
    384         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    385           feasible = false;
    386       } else if(rowActivity[i]>rowUpper[i]) {
    387         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    388           feasible = false;
    389       }
    390     }
    391     for (int i=0; i<numberIntegers; i++) {
    392       int iColumn = integerVariable[i];
    393       double value=newSolution[iColumn];
    394       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    395         feasible = false;
    396         break;
    397       }
    398     }
    399     if (feasible) {
    400       // new solution
    401       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    402       solutionValue = newSolutionValue;
    403       //printf("** Solution of %g found by CbcHeuristicDiveFractional\n",newSolutionValue);
    404       returnCode=1;
    405     } else {
    406       // Can easily happen
    407       //printf("Debug CbcHeuristicDiveFractional giving bad solution\n");
    408     }
    409   }
    410 
    411   delete [] newSolution;
    412   delete [] columnFixed;
    413   delete [] originalBound;
    414   delete [] fixedAtLowerBound;
    415   delete [] rowActivity;
    416   delete solver;
    417   return returnCode;
    418 }
    419 
    420 // update model
    421 void CbcHeuristicDiveFractional::setModel(CbcModel * model)
    422 {
    423   model_ = model;
    424   // Get a copy of original matrix (and by row for rounding);
    425   assert(model_->solver());
    426   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    427   if (matrix) {
    428     matrix_ = *matrix;
    429     // make sure model okay for heuristic
    430     validate();
    431102  }
    432103}
    433 
    434 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    435 void
    436 CbcHeuristicDiveFractional::validate()
    437 {
    438   if (model_&&when()<10) {
    439     if (model_->numberIntegers()!=
    440         model_->numberObjects())
    441       setWhen(0);
    442   }
    443 
    444   int numberIntegers = model_->numberIntegers();
    445   const int * integerVariable = model_->integerVariable();
    446   delete [] downLocks_;
    447   delete [] upLocks_;
    448   downLocks_ = new unsigned short [numberIntegers];
    449   upLocks_ = new unsigned short [numberIntegers];
    450   // Column copy
    451   const double * element = matrix_.getElements();
    452   const int * row = matrix_.getIndices();
    453   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    454   const int * columnLength = matrix_.getVectorLengths();
    455   const double * rowLower = model_->solver()->getRowLower();
    456   const double * rowUpper = model_->solver()->getRowUpper();
    457   for (int i=0;i<numberIntegers;i++) {
    458     int iColumn = integerVariable[i];
    459     int down=0;
    460     int up=0;
    461     if (columnLength[iColumn]>65535) {
    462       setWhen(0);
    463       break; // unlikely to work
    464     }
    465     for (CoinBigIndex j=columnStart[iColumn];
    466          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    467       int iRow=row[j];
    468       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    469         up++;
    470         down++;
    471       } else if (element[j]>0.0) {
    472         if (rowUpper[iRow]<1.0e20)
    473           up++;
    474         else
    475           down++;
    476       } else {
    477         if (rowLower[iRow]>-1.0e20)
    478           up++;
    479         else
    480           down++;
    481       }
    482     }
    483     downLocks_[i] = (unsigned short) down;
    484     upLocks_[i] = (unsigned short) up;
    485   }
    486 }
  • trunk/Cbc/src/CbcHeuristicDiveFractional.hpp

    r868 r912  
    44#define CbcHeuristicDiveFractional_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveFractional class
    99 */
    1010
    11 class CbcHeuristicDiveFractional : public CbcHeuristic {
     11class CbcHeuristicDiveFractional : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
    37 
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does Fractional Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     35  /// Selects the next variable to branch on
     36  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     37                                      const double* newSolution,
     38                                      int& bestColumn,
     39                                      int& bestRound);
    8640
    8741};
  • trunk/Cbc/src/CbcHeuristicDiveGuided.cpp

    r871 r912  
    88#include "CbcHeuristicDiveGuided.hpp"
    99#include "CbcStrategy.hpp"
    10 #include  "CoinTime.hpp"
    1110
    1211// Default Constructor
    1312CbcHeuristicDiveGuided::CbcHeuristicDiveGuided()
    14   :CbcHeuristic()
     13  :CbcHeuristicDive()
    1514{
    16   // matrix and row copy will automatically be empty
    17   downLocks_ =NULL;
    18   upLocks_ = NULL;
    19   percentageToFix_ = 0.2;
    20   maxIterations_ = 100;
    21   maxTime_ = 60;
    2215}
    2316
    2417// Constructor from model
    2518CbcHeuristicDiveGuided::CbcHeuristicDiveGuided(CbcModel & model)
    26   :CbcHeuristic(model)
     19  :CbcHeuristicDive(model)
    2720{
    28   downLocks_ =NULL;
    29   upLocks_ = NULL;
    30   // Get a copy of original matrix
    31   assert(model.solver());
    32   // model may have empty matrix - wait until setModel
    33   const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
    34   if (matrix) {
    35     matrix_ = *matrix;
    36   }
    37   percentageToFix_ = 0.2;
    38   maxIterations_ = 100;
    39   maxTime_ = 60;
    4021}
    4122
     
    4324CbcHeuristicDiveGuided::~CbcHeuristicDiveGuided ()
    4425{
    45   delete [] downLocks_;
    46   delete [] upLocks_;
    4726}
    4827
     
    6241  fprintf(fp,"3  CbcHeuristicDiveGuided heuristicDiveGuided(*cbcModel);\n");
    6342  CbcHeuristic::generateCpp(fp,"heuristicDiveGuided");
    64   if (percentageToFix_!=other.percentageToFix_)
    65     fprintf(fp,"3  heuristicDiveGuided.setPercentageToFix(%.2f);\n",percentageToFix_);
    66   else
    67     fprintf(fp,"4  heuristicDiveGuided.setPercentageToFix(%.2f);\n",percentageToFix_);
    68   if (maxIterations_!=other.maxIterations_)
    69     fprintf(fp,"3  heuristicDiveGuided.setMaxIterations(%d);\n",maxIterations_);
    70   else
    71     fprintf(fp,"4  heuristicDiveGuided.setMaxIterations(%d);\n",maxIterations_);
    72   if (maxTime_!=other.maxTime_)
    73     fprintf(fp,"3  heuristicDiveGuided.setMaxTime(%.2f);\n",maxTime_);
    74   else
    75     fprintf(fp,"4  heuristicDiveGuided.setMaxTime(%.2f);\n",maxTime_);
    7643  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveGuided);\n");
    7744}
     
    8047CbcHeuristicDiveGuided::CbcHeuristicDiveGuided(const CbcHeuristicDiveGuided & rhs)
    8148:
    82   CbcHeuristic(rhs),
    83   matrix_(rhs.matrix_),
    84   percentageToFix_(rhs.percentageToFix_),
    85   maxIterations_(rhs.maxIterations_),
    86   maxTime_(rhs.maxTime_)
     49  CbcHeuristicDive(rhs)
    8750{
    88   int numberIntegers = model_->numberIntegers();
    89   if (rhs.downLocks_) {
    90     int numberIntegers = model_->numberIntegers();
    91     downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    92     upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    93   } else {
    94     downLocks_ = NULL;
    95     upLocks_ = NULL;
    96   }
    9751}
    9852
     
    10256{
    10357  if (this!=&rhs) {
    104     CbcHeuristic::operator=(rhs);
    105     matrix_ = rhs.matrix_;
    106     percentageToFix_ = rhs.percentageToFix_;
    107     maxIterations_ = rhs.maxIterations_;
    108     maxTime_ = rhs.maxTime_;
    109     delete [] downLocks_;
    110     delete [] upLocks_;
    111     if (rhs.downLocks_) {
    112       int numberIntegers = model_->numberIntegers();
    113       downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    114       upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    115     } else {
    116       downLocks_ = NULL;
    117       upLocks_ = NULL;
    118     }
     58    CbcHeuristicDive::operator=(rhs);
    11959  }
    12060  return *this;
    12161}
    12262
    123 // Resets stuff if model changes
    124 void
    125 CbcHeuristicDiveGuided::resetModel(CbcModel * model)
     63bool
     64CbcHeuristicDiveGuided::canHeuristicRun()
    12665{
    127   model_=model;
    128   // Get a copy of original matrix
    129   assert(model_->solver());
    130   // model may have empty matrix - wait until setModel
    131   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    132   if (matrix) {
    133     matrix_ = *matrix;
    134     validate();
    135   }
     66  double* bestIntegerSolution = model_->bestSolution();
     67  if(bestIntegerSolution == NULL)
     68    return false; // no integer solution available. Switch off heuristic
     69
     70  return CbcHeuristicDive::canHeuristicRun();
    13671}
    13772
    138 // See if dive guided will give better solution
    139 // Sets value of solution
    140 // Returns 1 if solution, 0 if not
    141 int
    142 CbcHeuristicDiveGuided::solution(double & solutionValue,
    143                                  double * betterSolution)
     73void
     74CbcHeuristicDiveGuided::selectVariableToBranch(OsiSolverInterface* solver,
     75                                                   const double* newSolution,
     76                                                   int& bestColumn,
     77                                                   int& bestRound)
    14478{
     79  double* bestIntegerSolution = model_->bestSolution();
    14580
    146   // See if to do
    147   if (!when()||(when()%10==1&&model_->phase()!=1)||
    148       (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
    149     return 0; // switched off
    150 
    151   double* bestIntegerSolution = model_->bestSolution();
    152   if(bestIntegerSolution == NULL)
    153     return 0; // no integer solution available. Switch off heuristic
    154 
    155   double time1 = CoinCpuTime();
    156 
    157   OsiSolverInterface * solver = model_->solver()->clone();
    158   const double * lower = solver->getColLower();
    159   const double * upper = solver->getColUpper();
    160   const double * rowLower = solver->getRowLower();
    161   const double * rowUpper = solver->getRowUpper();
    162   const double * solution = solver->getColSolution();
    163   const double * objective = solver->getObjCoefficients();
    164   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    165   double primalTolerance;
    166   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    167 
    168   int numberRows = matrix_.getNumRows();
    169   assert (numberRows<=solver->getNumRows());
    17081  int numberIntegers = model_->numberIntegers();
    17182  const int * integerVariable = model_->integerVariable();
    172   double direction = solver->getObjSense(); // 1 for min, -1 for max
    173   double newSolutionValue = direction*solver->getObjValue();
    174   int returnCode = 0;
    175   // Column copy
    176   const double * element = matrix_.getElements();
    177   const int * row = matrix_.getIndices();
    178   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    179   const int * columnLength = matrix_.getVectorLengths();
     83  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    18084
    181   // Get solution array for heuristic solution
    182   int numberColumns = solver->getNumCols();
    183   double * newSolution = new double [numberColumns];
    184   memcpy(newSolution,solution,numberColumns*sizeof(double));
    185 
    186   // vectors to store the latest variables fixed at their bounds
    187   int* columnFixed = new int [numberIntegers];
    188   double* originalBound = new double [numberIntegers];
    189   bool * fixedAtLowerBound = new bool [numberIntegers];
    190 
    191   const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
    192 
    193   // count how many fractional variables
    194   int numberFractionalVariables = 0;
     85  bestColumn = -1;
     86  bestRound = -1; // -1 rounds down, +1 rounds up
     87  double bestFraction = DBL_MAX;
    19588  for (int i=0; i<numberIntegers; i++) {
    19689    int iColumn = integerVariable[i];
    19790    double value=newSolution[iColumn];
     91    double fraction=value-floor(value);
     92    int round = 0;
    19893    if (fabs(floor(value+0.5)-value)>integerTolerance) {
    199       numberFractionalVariables++;
    200     }
    201   }
     94      if(downLocks_[i]>0&&upLocks_[i]>0) {
     95        // the variable cannot be rounded
     96        if(value >= bestIntegerSolution[iColumn])
     97          round = -1;
     98        else {
     99          round = 1;
     100          fraction = 1.0 - fraction;
     101        }
     102       
     103        // if variable is not binary, penalize it
     104        if(!solver->isBinary(iColumn))
     105          fraction *= 1000.0;
    202106
    203   int iteration = 0;
    204   while(numberFractionalVariables) {
    205     iteration++;
    206 
    207     // select a fractional variable to bound
    208     double bestFraction = DBL_MAX;
    209     int bestColumn = -1;
    210     int bestRound = -1; // -1 rounds down, +1 rounds up
    211     int numberAtBoundFixed = 0;
    212     bool canRoundSolution = true;
    213     for (int i=0; i<numberIntegers; i++) {
    214       int iColumn = integerVariable[i];
    215       double value=newSolution[iColumn];
    216       double fraction=value-floor(value);
    217       int round = 0;
    218       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    219         if(downLocks_[i]>0&&upLocks_[i]>0) {
    220           // the variable cannot be rounded
    221           canRoundSolution = false;
    222           if(value >= bestIntegerSolution[iColumn])
    223             round = -1;
    224           else {
    225             round = 1;
    226             fraction = 1.0 - fraction;
    227           }
    228 
    229           // if variable is not binary, penalize it
    230           if(!solver->isBinary(iColumn))
    231             fraction *= 1000.0;
    232 
    233           if(fraction < bestFraction) {
    234             bestColumn = iColumn;
    235             bestFraction = fraction;
    236             bestRound = round;
    237           }
    238         }
    239       }
    240       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    241         // fix the variable at one of its bounds
    242         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    243             lower[iColumn] != upper[iColumn]) {
    244           columnFixed[numberAtBoundFixed] = iColumn;
    245           originalBound[numberAtBoundFixed] = upper[iColumn];
    246           fixedAtLowerBound[numberAtBoundFixed] = true;
    247           solver->setColUpper(iColumn, lower[iColumn]);
    248           numberAtBoundFixed++;
    249         }
    250         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    251             lower[iColumn] != upper[iColumn]) {
    252           columnFixed[numberAtBoundFixed] = iColumn;
    253           originalBound[numberAtBoundFixed] = lower[iColumn];
    254           fixedAtLowerBound[numberAtBoundFixed] = false;
    255           solver->setColLower(iColumn, upper[iColumn]);
    256           numberAtBoundFixed++;
     107        if(fraction < bestFraction) {
     108          bestColumn = iColumn;
     109          bestFraction = fraction;
     110          bestRound = round;
    257111        }
    258112      }
    259113    }
    260 
    261     if(canRoundSolution) {
    262       // Round all the fractional variables
    263       for (int i=0; i<numberIntegers; i++) {
    264         int iColumn = integerVariable[i];
    265         double value=newSolution[iColumn];
    266         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    267           if(downLocks_[i]==0 || upLocks_[i]==0) {
    268             if(downLocks_[i]==0 && upLocks_[i]==0) {
    269               if(direction * objective[iColumn] >= 0.0)
    270                 newSolution[iColumn] = floor(value);
    271               else
    272                 newSolution[iColumn] = ceil(value);
    273             }
    274             else if(downLocks_[i]==0)
    275               newSolution[iColumn] = floor(value);
    276             else
    277               newSolution[iColumn] = ceil(value);
    278           }
    279           else
    280             break;
    281         }
    282       }
    283       break;
    284     }
    285 
    286 
    287     double originalBoundBestColumn;
    288     if(bestColumn >= 0) {
    289       if(bestRound < 0) {
    290         originalBoundBestColumn = upper[bestColumn];
    291         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    292       }
    293       else {
    294         originalBoundBestColumn = lower[bestColumn];
    295         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    296       }
    297     } else {
    298       break;
    299     }
    300     int originalBestRound = bestRound;
    301     while (1) {
    302 
    303       solver->resolve();
    304 
    305       if(!solver->isProvenOptimal()) {
    306         if(numberAtBoundFixed > 0) {
    307           // Remove the bound fix for variables that were at bounds
    308           for(int i=0; i<numberAtBoundFixed; i++) {
    309             int iColFixed = columnFixed[i];
    310             if(fixedAtLowerBound[i])
    311               solver->setColUpper(iColFixed, originalBound[i]);
    312             else
    313               solver->setColLower(iColFixed, originalBound[i]);
    314           }
    315           numberAtBoundFixed = 0;
    316         }
    317         else if(bestRound == originalBestRound) {
    318           bestRound *= (-1);
    319           if(bestRound < 0) {
    320             solver->setColLower(bestColumn, originalBoundBestColumn);
    321             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    322           }
    323           else {
    324             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    325             solver->setColUpper(bestColumn, originalBoundBestColumn);
    326           }
    327         }
    328         else
    329           break;
    330       }
    331       else
    332         break;
    333     }
    334 
    335     if(!solver->isProvenOptimal())
    336       break;
    337 
    338     if(iteration > maxIterations_) {
    339       break;
    340     }
    341 
    342     if(CoinCpuTime()-time1 > maxTime_) {
    343       break;
    344     }
    345 
    346     memcpy(newSolution,solution,numberColumns*sizeof(double));
    347     numberFractionalVariables = 0;
    348     for (int i=0; i<numberIntegers; i++) {
    349       int iColumn = integerVariable[i];
    350       double value=newSolution[iColumn];
    351       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    352         numberFractionalVariables++;
    353       }
    354     }
    355 
    356   }
    357 
    358 
    359   double * rowActivity = new double[numberRows];
    360   memset(rowActivity,0,numberRows*sizeof(double));
    361 
    362   // re-compute new solution value
    363   double objOffset=0.0;
    364   solver->getDblParam(OsiObjOffset,objOffset);
    365   newSolutionValue = -objOffset;
    366   for (int i=0 ; i<numberColumns ; i++ )
    367     newSolutionValue += objective[i]*newSolution[i];
    368   newSolutionValue *= direction;
    369     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    370   if (newSolutionValue<solutionValue) {
    371     // paranoid check
    372     memset(rowActivity,0,numberRows*sizeof(double));
    373     for (int i=0;i<numberColumns;i++) {
    374       int j;
    375       double value = newSolution[i];
    376       if (value) {
    377         for (j=columnStart[i];
    378              j<columnStart[i]+columnLength[i];j++) {
    379           int iRow=row[j];
    380           rowActivity[iRow] += value*element[j];
    381         }
    382       }
    383     }
    384     // check was approximately feasible
    385     bool feasible=true;
    386     for (int i=0;i<numberRows;i++) {
    387       if(rowActivity[i]<rowLower[i]) {
    388         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    389           feasible = false;
    390       } else if(rowActivity[i]>rowUpper[i]) {
    391         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    392           feasible = false;
    393       }
    394     }
    395     for (int i=0; i<numberIntegers; i++) {
    396       int iColumn = integerVariable[i];
    397       double value=newSolution[iColumn];
    398       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    399         feasible = false;
    400         break;
    401       }
    402     }
    403     if (feasible) {
    404       // new solution
    405       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    406       solutionValue = newSolutionValue;
    407       //printf("** Solution of %g found by CbcHeuristicDiveGuided\n",newSolutionValue);
    408       returnCode=1;
    409     } else {
    410       // Can easily happen
    411       //printf("Debug CbcHeuristicDiveGuided giving bad solution\n");
    412     }
    413   }
    414 
    415   delete [] newSolution;
    416   delete [] columnFixed;
    417   delete [] originalBound;
    418   delete [] fixedAtLowerBound;
    419   delete [] rowActivity;
    420   delete solver;
    421   return returnCode;
    422 }
    423 
    424 // update model
    425 void CbcHeuristicDiveGuided::setModel(CbcModel * model)
    426 {
    427   model_ = model;
    428   // Get a copy of original matrix (and by row for rounding);
    429   assert(model_->solver());
    430   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    431   if (matrix) {
    432     matrix_ = *matrix;
    433     // make sure model okay for heuristic
    434     validate();
    435114  }
    436115}
    437 
    438 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    439 void
    440 CbcHeuristicDiveGuided::validate()
    441 {
    442   if (model_&&when()<10) {
    443     if (model_->numberIntegers()!=
    444         model_->numberObjects())
    445       setWhen(0);
    446   }
    447 
    448   int numberIntegers = model_->numberIntegers();
    449   const int * integerVariable = model_->integerVariable();
    450   delete [] downLocks_;
    451   delete [] upLocks_;
    452   downLocks_ = new unsigned short [numberIntegers];
    453   upLocks_ = new unsigned short [numberIntegers];
    454   // Column copy
    455   const double * element = matrix_.getElements();
    456   const int * row = matrix_.getIndices();
    457   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    458   const int * columnLength = matrix_.getVectorLengths();
    459   const double * rowLower = model_->solver()->getRowLower();
    460   const double * rowUpper = model_->solver()->getRowUpper();
    461   for (int i=0;i<numberIntegers;i++) {
    462     int iColumn = integerVariable[i];
    463     int down=0;
    464     int up=0;
    465     if (columnLength[iColumn]>65535) {
    466       setWhen(0);
    467       break; // unlikely to work
    468     }
    469     for (CoinBigIndex j=columnStart[iColumn];
    470          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    471       int iRow=row[j];
    472       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    473         up++;
    474         down++;
    475       } else if (element[j]>0.0) {
    476         if (rowUpper[iRow]<1.0e20)
    477           up++;
    478         else
    479           down++;
    480       } else {
    481         if (rowLower[iRow]>-1.0e20)
    482           up++;
    483         else
    484           down++;
    485       }
    486     }
    487     downLocks_[i] = (unsigned short) down;
    488     upLocks_[i] = (unsigned short) up;
    489   }
    490 }
  • trunk/Cbc/src/CbcHeuristicDiveGuided.hpp

    r868 r912  
    44#define CbcHeuristicDiveGuided_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveGuided class
    99 */
    1010
    11 class CbcHeuristicDiveGuided : public CbcHeuristic {
     11class CbcHeuristicDiveGuided : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
     35  /// Tests if the heuristic can run
     36  virtual bool canHeuristicRun();
    3737
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does Guided Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     38  /// Selects the next variable to branch on
     39  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     40                                      const double* newSolution,
     41                                      int& bestColumn,
     42                                      int& bestRound);
    8643
    8744};
  • trunk/Cbc/src/CbcHeuristicDiveVectorLength.cpp

    r871 r912  
    88#include "CbcHeuristicDiveVectorLength.hpp"
    99#include "CbcStrategy.hpp"
    10 #include  "CoinTime.hpp"
    1110
    1211// Default Constructor
    1312CbcHeuristicDiveVectorLength::CbcHeuristicDiveVectorLength()
    14   :CbcHeuristic()
     13  :CbcHeuristicDive()
    1514{
    16   // matrix and row copy will automatically be empty
    17   downLocks_ =NULL;
    18   upLocks_ = NULL;
    19   percentageToFix_ = 0.2;
    20   maxIterations_ = 100;
    21   maxTime_ = 60;
    2215}
    2316
    2417// Constructor from model
    2518CbcHeuristicDiveVectorLength::CbcHeuristicDiveVectorLength(CbcModel & model)
    26   :CbcHeuristic(model)
     19  :CbcHeuristicDive(model)
    2720{
    28   // Get a copy of original matrix
    29   assert(model.solver());
    30   downLocks_ =NULL;
    31   upLocks_ = NULL;
    32   // model may have empty matrix - wait until setModel
    33   const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
    34   if (matrix) {
    35     matrix_ = *matrix;
    36   }
    37   percentageToFix_ = 0.2;
    38   maxIterations_ = 100;
    39   maxTime_ = 60;
    4021}
    4122
     
    4324CbcHeuristicDiveVectorLength::~CbcHeuristicDiveVectorLength ()
    4425{
    45   delete [] downLocks_;
    46   delete [] upLocks_;
    4726}
    4827
     
    6241  fprintf(fp,"3  CbcHeuristicDiveVectorLength heuristicDiveVectorLength(*cbcModel);\n");
    6342  CbcHeuristic::generateCpp(fp,"heuristicDiveVectorLength");
    64   if (percentageToFix_!=other.percentageToFix_)
    65     fprintf(fp,"3  heuristicDiveVectorLength.setPercentageToFix(%.2f);\n",percentageToFix_);
    66   else
    67     fprintf(fp,"4  heuristicDiveVectorLength.setPercentageToFix(%.2f);\n",percentageToFix_);
    68   if (maxIterations_!=other.maxIterations_)
    69     fprintf(fp,"3  heuristicDiveVectorLength.setMaxIterations(%d);\n",maxIterations_);
    70   else
    71     fprintf(fp,"4  heuristicDiveVectorLength.setMaxIterations(%d);\n",maxIterations_);
    72   if (maxTime_!=other.maxTime_)
    73     fprintf(fp,"3  heuristicDiveVectorLength.setMaxTime(%.2f);\n",maxTime_);
    74   else
    75     fprintf(fp,"4  heuristicDiveVectorLength.setMaxTime(%.2f);\n",maxTime_);
    7643  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveVectorLength);\n");
    7744}
     
    8047CbcHeuristicDiveVectorLength::CbcHeuristicDiveVectorLength(const CbcHeuristicDiveVectorLength & rhs)
    8148:
    82   CbcHeuristic(rhs),
    83   matrix_(rhs.matrix_),
    84   percentageToFix_(rhs.percentageToFix_),
    85   maxIterations_(rhs.maxIterations_),
    86   maxTime_(rhs.maxTime_)
     49  CbcHeuristicDive(rhs)
    8750{
    88   if (rhs.downLocks_) {
    89     int numberIntegers = model_->numberIntegers();
    90     downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    91     upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    92   } else {
    93     downLocks_ = NULL;
    94     upLocks_ = NULL;
    95   }
    9651}
    9752
     
    10156{
    10257  if (this!=&rhs) {
    103     CbcHeuristic::operator=(rhs);
    104     matrix_ = rhs.matrix_;
    105     percentageToFix_ = rhs.percentageToFix_;
    106     maxIterations_ = rhs.maxIterations_;
    107     maxTime_ = rhs.maxTime_;
    108     delete [] downLocks_;
    109     delete [] upLocks_;
    110     int numberIntegers = model_->numberIntegers();
    111     if (rhs.downLocks_) {
    112       downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    113       upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    114     } else {
    115       downLocks_ = NULL;
    116       upLocks_ = NULL;
    117     }
     58    CbcHeuristicDive::operator=(rhs);
    11859  }
    11960  return *this;
    12061}
    12162
    122 // Resets stuff if model changes
    123 void
    124 CbcHeuristicDiveVectorLength::resetModel(CbcModel * model)
     63void
     64CbcHeuristicDiveVectorLength::selectVariableToBranch(OsiSolverInterface* solver,
     65                                                   const double* newSolution,
     66                                                   int& bestColumn,
     67                                                   int& bestRound)
    12568{
    126   model_=model;
    127   // Get a copy of original matrix
    128   assert(model_->solver());
    129   // model may have empty matrix - wait until setModel
    130   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    131   if (matrix) {
    132     matrix_ = *matrix;
    133     validate();
    134   }
    135 }
     69  const double * objective = solver->getObjCoefficients();
     70  double direction = solver->getObjSense(); // 1 for min, -1 for max
    13671
    137 // See if dive fractional will give better solution
    138 // Sets value of solution
    139 // Returns 1 if solution, 0 if not
    140 int
    141 CbcHeuristicDiveVectorLength::solution(double & solutionValue,
    142                                      double * betterSolution)
    143 {
    144 
    145   // See if to do
    146   if (!when()||(when()%10==1&&model_->phase()!=1)||
    147       (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
    148     return 0; // switched off
    149 
    150   double time1 = CoinCpuTime();
    151 
    152   OsiSolverInterface * solver = model_->solver()->clone();
    153   const double * lower = solver->getColLower();
    154   const double * upper = solver->getColUpper();
    155   const double * rowLower = solver->getRowLower();
    156   const double * rowUpper = solver->getRowUpper();
    157   const double * solution = solver->getColSolution();
    158   const double * objective = solver->getObjCoefficients();
    159   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    160   double primalTolerance;
    161   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    162 
    163   int numberRows = matrix_.getNumRows();
    164   assert (numberRows<=solver->getNumRows());
     72  const int * columnLength = matrix_.getVectorLengths();
    16573  int numberIntegers = model_->numberIntegers();
    16674  const int * integerVariable = model_->integerVariable();
    167   double direction = solver->getObjSense(); // 1 for min, -1 for max
    168   double newSolutionValue = direction*solver->getObjValue();
    169   int returnCode = 0;
    170   // Column copy
    171   const double * element = matrix_.getElements();
    172   const int * row = matrix_.getIndices();
    173   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    174   const int * columnLength = matrix_.getVectorLengths();
     75  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    17576
    176   // Get solution array for heuristic solution
    177   int numberColumns = solver->getNumCols();
    178   double * newSolution = new double [numberColumns];
    179   memcpy(newSolution,solution,numberColumns*sizeof(double));
    180 
    181   // vectors to store the latest variables fixed at their bounds
    182   int* columnFixed = new int [numberIntegers];
    183   double* originalBound = new double [numberIntegers];
    184   bool * fixedAtLowerBound = new bool [numberIntegers];
    185 
    186   const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
    187 
    188   // count how many fractional variables
    189   int numberFractionalVariables = 0;
     77  bestColumn = -1;
     78  bestRound = -1; // -1 rounds down, +1 rounds up
     79  double bestScore = DBL_MAX;
    19080  for (int i=0; i<numberIntegers; i++) {
    19181    int iColumn = integerVariable[i];
    19282    double value=newSolution[iColumn];
     83    double fraction=value-floor(value);
     84    int round = 0;
    19385    if (fabs(floor(value+0.5)-value)>integerTolerance) {
    194       numberFractionalVariables++;
    195     }
    196   }
     86      if(downLocks_[i]>0&&upLocks_[i]>0) {
     87        // the variable cannot be rounded
     88        double obj = direction * objective[iColumn];
     89        if(obj >= 0.0)
     90          round = 1; // round up
     91        else
     92          round = -1; // round down
     93        double objDelta;
     94        if(round == 1)
     95          objDelta = (1.0 - fraction) * obj;
     96        else
     97          objDelta = - fraction * obj;
     98       
     99        // we want the smaller score
     100        double score = objDelta / ((double) columnLength[iColumn] + 1.0);
    197101
    198   int iteration = 0;
    199   while(numberFractionalVariables) {
    200     iteration++;
     102        // if variable is not binary, penalize it
     103        if(!solver->isBinary(iColumn))
     104          score *= 1000.0;
    201105
    202     // select a fractional variable to bound
    203     int bestColumn = -1;
    204     int bestRound = -1; // -1 rounds down, +1 rounds up
    205     double bestScore = DBL_MAX;
    206     int numberAtBoundFixed = 0;
    207     bool canRoundSolution = true;
    208     for (int i=0; i<numberIntegers; i++) {
    209       int iColumn = integerVariable[i];
    210       double value=newSolution[iColumn];
    211       double fraction=value-floor(value);
    212       int round = 0;
    213       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    214         if(downLocks_[i]>0 && upLocks_[i]>0) {
    215           // the variable cannot be rounded
    216           canRoundSolution = false;
    217           double obj = direction * objective[iColumn];
    218           if(obj >= 0.0)
    219             round = 1; // round up
    220           else
    221             round = -1; // round down
    222           double objDelta;
    223           if(round == 1)
    224             objDelta = (1.0 - fraction) * obj;
    225           else
    226             objDelta = - fraction * obj;
    227 
    228           // we want the smaller score
    229           double score = objDelta / ((double) columnLength[iColumn] + 1.0);
    230 
    231           // if variable is not binary, penalize it
    232           if(!solver->isBinary(iColumn))
    233             score *= 1000.0;
    234 
    235           if(score < bestScore) {
    236             bestColumn = iColumn;
    237             bestScore = score;
    238             bestRound = round;
    239           }
    240         }
    241       }
    242       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    243         // fix the variable at one of its bounds
    244         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    245             lower[iColumn] != upper[iColumn]) {
    246           columnFixed[numberAtBoundFixed] = iColumn;
    247           originalBound[numberAtBoundFixed] = upper[iColumn];
    248           fixedAtLowerBound[numberAtBoundFixed] = true;
    249           solver->setColUpper(iColumn, lower[iColumn]);
    250           numberAtBoundFixed++;
    251         }
    252         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    253             lower[iColumn] != upper[iColumn]) {
    254           columnFixed[numberAtBoundFixed] = iColumn;
    255           originalBound[numberAtBoundFixed] = lower[iColumn];
    256           fixedAtLowerBound[numberAtBoundFixed] = false;
    257           solver->setColLower(iColumn, upper[iColumn]);
    258           numberAtBoundFixed++;
     106        if(score < bestScore) {
     107          bestColumn = iColumn;
     108          bestScore = score;
     109          bestRound = round;
    259110        }
    260111      }
    261112    }
    262 
    263     if(canRoundSolution) {
    264       // Round all the fractional variables
    265       for (int i=0; i<numberIntegers; i++) {
    266         int iColumn = integerVariable[i];
    267         double value=newSolution[iColumn];
    268         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    269           if(downLocks_[i]==0 || upLocks_[i]==0) {
    270             if(downLocks_[i]==0 && upLocks_[i]==0) {
    271               if(direction * objective[iColumn] >= 0.0)
    272                 newSolution[iColumn] = floor(value);
    273               else
    274                 newSolution[iColumn] = ceil(value);
    275             }
    276             else if(downLocks_[i]==0)
    277               newSolution[iColumn] = floor(value);
    278             else
    279               newSolution[iColumn] = ceil(value);
    280           }
    281           else
    282             break;
    283         }
    284       }
    285       break;
    286     }
    287 
    288 
    289     double originalBoundBestColumn;
    290     if(bestColumn >= 0) {
    291       if(bestRound < 0) {
    292         originalBoundBestColumn = upper[bestColumn];
    293         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    294       }
    295       else {
    296         originalBoundBestColumn = lower[bestColumn];
    297         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    298       }
    299     } else {
    300       break;
    301     }
    302     int originalBestRound = bestRound;
    303     while (1) {
    304 
    305       solver->resolve();
    306 
    307       if(!solver->isProvenOptimal()) {
    308         if(numberAtBoundFixed > 0) {
    309           // Remove the bound fix for variables that were at bounds
    310           for(int i=0; i<numberAtBoundFixed; i++) {
    311             int iColFixed = columnFixed[i];
    312             if(fixedAtLowerBound[i])
    313               solver->setColUpper(iColFixed, originalBound[i]);
    314             else
    315               solver->setColLower(iColFixed, originalBound[i]);
    316           }
    317           numberAtBoundFixed = 0;
    318         }
    319         else if(bestRound == originalBestRound) {
    320           bestRound *= (-1);
    321           if(bestRound < 0) {
    322             solver->setColLower(bestColumn, originalBoundBestColumn);
    323             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    324           }
    325           else {
    326             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    327             solver->setColUpper(bestColumn, originalBoundBestColumn);
    328           }
    329         }
    330         else
    331           break;
    332       }
    333       else
    334         break;
    335     }
    336 
    337     if(!solver->isProvenOptimal())
    338       break;
    339 
    340     if(iteration > maxIterations_) {
    341       break;
    342     }
    343 
    344     if(CoinCpuTime()-time1 > maxTime_) {
    345       break;
    346     }
    347 
    348     memcpy(newSolution,solution,numberColumns*sizeof(double));
    349     numberFractionalVariables = 0;
    350     for (int i=0; i<numberIntegers; i++) {
    351       int iColumn = integerVariable[i];
    352       double value=newSolution[iColumn];
    353       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    354         numberFractionalVariables++;
    355       }
    356     }
    357 
    358   }
    359 
    360 
    361   double * rowActivity = new double[numberRows];
    362   memset(rowActivity,0,numberRows*sizeof(double));
    363 
    364   // re-compute new solution value
    365   double objOffset=0.0;
    366   solver->getDblParam(OsiObjOffset,objOffset);
    367   newSolutionValue = -objOffset;
    368   for (int i=0 ; i<numberColumns ; i++ )
    369     newSolutionValue += objective[i]*newSolution[i];
    370   newSolutionValue *= direction;
    371     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    372   if (newSolutionValue<solutionValue) {
    373     // paranoid check
    374     memset(rowActivity,0,numberRows*sizeof(double));
    375     for (int i=0;i<numberColumns;i++) {
    376       int j;
    377       double value = newSolution[i];
    378       if (value) {
    379         for (j=columnStart[i];
    380              j<columnStart[i]+columnLength[i];j++) {
    381           int iRow=row[j];
    382           rowActivity[iRow] += value*element[j];
    383         }
    384       }
    385     }
    386     // check was approximately feasible
    387     bool feasible=true;
    388     for (int i=0;i<numberRows;i++) {
    389       if(rowActivity[i]<rowLower[i]) {
    390         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    391           feasible = false;
    392       } else if(rowActivity[i]>rowUpper[i]) {
    393         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    394           feasible = false;
    395       }
    396     }
    397     for (int i=0; i<numberIntegers; i++) {
    398       int iColumn = integerVariable[i];
    399       double value=newSolution[iColumn];
    400       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    401         feasible = false;
    402         break;
    403       }
    404     }
    405     if (feasible) {
    406       // new solution
    407       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    408       solutionValue = newSolutionValue;
    409       //printf("** Solution of %g found by CbcHeuristicDiveVectorLength\n",newSolutionValue);
    410       returnCode=1;
    411     } else {
    412       // Can easily happen
    413       //printf("Debug CbcHeuristicDiveVectorLength giving bad solution\n");
    414     }
    415   }
    416 
    417   delete [] newSolution;
    418   delete [] columnFixed;
    419   delete [] originalBound;
    420   delete [] fixedAtLowerBound;
    421   delete [] rowActivity;
    422   delete solver;
    423   return returnCode;
    424 }
    425 
    426 // update model
    427 void CbcHeuristicDiveVectorLength::setModel(CbcModel * model)
    428 {
    429   model_ = model;
    430   // Get a copy of original matrix (and by row for rounding);
    431   assert(model_->solver());
    432   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    433   if (matrix) {
    434     matrix_ = *matrix;
    435     // make sure model okay for heuristic
    436     validate();
    437113  }
    438114}
    439 
    440 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    441 void
    442 CbcHeuristicDiveVectorLength::validate()
    443 {
    444   if (model_&&when()<10) {
    445     if (model_->numberIntegers()!=
    446         model_->numberObjects())
    447       setWhen(0);
    448   }
    449 
    450   int numberIntegers = model_->numberIntegers();
    451   const int * integerVariable = model_->integerVariable();
    452   delete [] downLocks_;
    453   delete [] upLocks_;
    454   downLocks_ = new unsigned short [numberIntegers];
    455   upLocks_ = new unsigned short [numberIntegers];
    456   // Column copy
    457   const double * element = matrix_.getElements();
    458   const int * row = matrix_.getIndices();
    459   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    460   const int * columnLength = matrix_.getVectorLengths();
    461   const double * rowLower = model_->solver()->getRowLower();
    462   const double * rowUpper = model_->solver()->getRowUpper();
    463   for (int i=0;i<numberIntegers;i++) {
    464     int iColumn = integerVariable[i];
    465     int down=0;
    466     int up=0;
    467     if (columnLength[iColumn]>65535) {
    468       setWhen(0);
    469       break; // unlikely to work
    470     }
    471     for (CoinBigIndex j=columnStart[iColumn];
    472          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    473       int iRow=row[j];
    474       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    475         up++;
    476         down++;
    477       } else if (element[j]>0.0) {
    478         if (rowUpper[iRow]<1.0e20)
    479           up++;
    480         else
    481           down++;
    482       } else {
    483         if (rowLower[iRow]>-1.0e20)
    484           up++;
    485         else
    486           down++;
    487       }
    488     }
    489     downLocks_[i] = (unsigned short) down;
    490     upLocks_[i] = (unsigned short) up;
    491   }
    492 }
  • trunk/Cbc/src/CbcHeuristicDiveVectorLength.hpp

    r868 r912  
    44#define CbcHeuristicDiveVectorLength_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveVectorLength class
    99 */
    1010
    11 class CbcHeuristicDiveVectorLength : public CbcHeuristic {
     11class CbcHeuristicDiveVectorLength : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
    37 
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does VectorLength Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     35  /// Selects the next variable to branch on
     36  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     37                                      const double* newSolution,
     38                                      int& bestColumn,
     39                                      int& bestRound);
    8640
    8741};
  • trunk/Cbc/src/CbcNode.cpp

    r881 r912  
    4646  numberPointingToThis_(0),
    4747  parent_(NULL),
     48  parentBranch_(NULL),
    4849  owner_(NULL),
    4950  numberCuts_(0),
     
    5859#endif
    5960}
     61
     62void
     63CbcNodeInfo::setParentBasedData()
     64{
     65  if (parent_) {
     66    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
     67    //parent_->increment();
     68    if (parent_->owner()) {
     69      const OsiBranchingObject* br = parent_->owner()->branchingObject();
     70      const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
     71      assert(cbcbr);
     72      parentBranch_ = cbcbr->clone();
     73      parentBranch_->previousBranch();
     74    }
     75  }
     76}
     77
     78#if 0
    6079// Constructor given parent
    6180CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent)
     
    6382  numberPointingToThis_(2),
    6483  parent_(parent),
     84  parentBranch_(NULL),
    6585  owner_(NULL),
    6686  numberCuts_(0),
     
    7494  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
    7595#endif
    76   if (parent_) {
    77     numberRows_ = parent_->numberRows_+parent_->numberCuts_;
    78     //parent_->increment();
    79   }
     96  setParentBasedData();
    8097}
     98#endif
     99
    81100// Copy Constructor
    82101CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs)
     
    84103  numberPointingToThis_(rhs.numberPointingToThis_),
    85104  parent_(rhs.parent_),
     105  parentBranch_(NULL),
    86106  owner_(rhs.owner_),
    87107  numberCuts_(rhs.numberCuts_),
     
    109129    numberCuts_=n;
    110130  }
     131  if (rhs.parentBranch_) {
     132    parentBranch_ = rhs.parentBranch_->clone();
     133  }
    111134}
    112135// Constructor given parent and owner
     
    115138  numberPointingToThis_(2),
    116139  parent_(parent),
     140  parentBranch_(NULL),
    117141  owner_(owner),
    118142  numberCuts_(0),
     
    126150  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
    127151#endif
    128   if (parent_) {
    129     numberRows_ = parent_->numberRows_+parent_->numberCuts_;
    130   }
     152  setParentBasedData();
    131153}
    132154
     
    161183    if (!numberLinks) delete parent_;
    162184  }
     185  delete parentBranch_;
    163186}
    164187
     
    479502CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model,
    480503                                 int numberRowsAtContinuous) :
    481   CbcNodeInfo()
     504  CbcNodeInfo(NULL, model->currentNode())
    482505{
    483506  OsiSolverInterface * solver = model->solver();
     
    654677CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs)
    655678
    656   : CbcNodeInfo(rhs.parent_)
     679  : CbcNodeInfo(rhs)
    657680
    658681{ basisDiff_ = rhs.basisDiff_->clone() ;
  • trunk/Cbc/src/CbcNode.hpp

    r854 r912  
    7474  CbcNodeInfo ( const CbcNodeInfo &);
    7575   
    76 
     76#if 0
    7777  /** Construct with parent
    7878
     
    8181  */
    8282  CbcNodeInfo (CbcNodeInfo * parent);
     83#endif
    8384   
    8485  /** Construct with parent and owner
     
    229230  inline void unmark()
    230231  { active_ &= ~8;}
     232
     233  /// Branching object for the parent
     234  inline const CbcBranchingObject * parentBranch() const
     235  { return parentBranch_;}
    231236protected:
    232237
     
    243248  CbcNodeInfo * parent_;
    244249
     250  /// Copy of the branching object of the parent when the node is created
     251  CbcBranchingObject * parentBranch_;
     252     
    245253  /// Owner
    246254  CbcNode * owner_;
     
    272280  */
    273281  int active_;
    274      
     282
    275283private:
    276284 
    277285  /// Illegal Assignment operator
    278286  CbcNodeInfo & operator=(const CbcNodeInfo& rhs);
    279  
     287
     288  /// routine common to constructors
     289  void setParentBasedData();
    280290};
    281291
     
    635645  inline void setBranchingObject(OsiBranchingObject * branchingObject)
    636646  { branch_ = branchingObject;}
    637   /// The node number
     647  /// The node number 
    638648  inline int nodeNumber() const
    639649  { return nodeNumber_;}
  • trunk/Cbc/src/Makefile.am

    r871 r912  
    3535        CbcFeasibilityBase.hpp \
    3636        CbcHeuristic.cpp CbcHeuristic.hpp \
     37        CbcHeuristicDive.cpp CbcHeuristicDive.hpp \
    3738        CbcHeuristicDiveCoefficient.cpp CbcHeuristicDiveCoefficient.hpp \
    3839        CbcHeuristicDiveFractional.cpp CbcHeuristicDiveFractional.hpp \
     
    307308        CbcFeasibilityBase.hpp \
    308309        CbcHeuristic.hpp \
     310        CbcHeuristicDive.hpp \
    309311        CbcHeuristicDiveCoefficient.hpp \
    310312        CbcHeuristicDiveFractional.hpp \
  • trunk/Cbc/src/Makefile.in

    r871 r912  
    162162        CbcCutGenerator.lo CbcEventHandler.lo CbcFathom.lo \
    163163        CbcFathomDynamicProgramming.lo CbcHeuristic.lo \
    164         CbcHeuristicDiveCoefficient.lo CbcHeuristicDiveFractional.lo \
    165         CbcHeuristicDiveGuided.lo CbcHeuristicDiveVectorLength.lo \
    166         CbcHeuristicFPump.lo CbcHeuristicGreedy.lo \
    167         CbcHeuristicLocal.lo CbcHeuristicRINS.lo CbcMessage.lo \
    168         CbcModel.lo CbcNode.lo CbcStatistics.lo CbcStrategy.lo \
    169         CbcTree.lo CbcTreeLocal.lo
     164        CbcHeuristicDive.lo CbcHeuristicDiveCoefficient.lo \
     165        CbcHeuristicDiveFractional.lo CbcHeuristicDiveGuided.lo \
     166        CbcHeuristicDiveVectorLength.lo CbcHeuristicFPump.lo \
     167        CbcHeuristicGreedy.lo CbcHeuristicLocal.lo CbcHeuristicRINS.lo \
     168        CbcMessage.lo CbcModel.lo CbcNode.lo CbcStatistics.lo \
     169        CbcStrategy.lo CbcTree.lo CbcTreeLocal.lo
    170170libCbc_la_OBJECTS = $(am_libCbc_la_OBJECTS)
    171171libCbcSolver_la_LIBADD =
     
    510510        CbcFeasibilityBase.hpp \
    511511        CbcHeuristic.cpp CbcHeuristic.hpp \
     512        CbcHeuristicDive.cpp CbcHeuristicDive.hpp \
    512513        CbcHeuristicDiveCoefficient.cpp CbcHeuristicDiveCoefficient.hpp \
    513514        CbcHeuristicDiveFractional.cpp CbcHeuristicDiveFractional.hpp \
     
    650651        CbcFeasibilityBase.hpp \
    651652        CbcHeuristic.hpp \
     653        CbcHeuristicDive.hpp \
    652654        CbcHeuristicDiveCoefficient.hpp \
    653655        CbcHeuristicDiveFractional.hpp \
     
    806808@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcGeneric.Po@am__quote@
    807809@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristic.Plo@am__quote@
     810@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristicDive.Plo@am__quote@
    808811@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristicDiveCoefficient.Plo@am__quote@
    809812@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristicDiveFractional.Plo@am__quote@
Note: See TracChangeset for help on using the changeset viewer.