Changeset 995 for branches


Ignore:
Timestamp:
Jun 29, 2008 6:49:32 AM (11 years ago)
Author:
ladanyi
Message:

seems to work

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/dynamicbranching/dynamicbranching.cpp

    r994 r995  
    7979  double* getColLower;
    8080  double* getColUpper;
     81  double* rowDualRay;
     82  double* colDualRay;
    8183  double* getColSolution; // FIXME Not needed, just for debugging
    82   double* getObjCoefficients;
    8384  double yb_plus_rl_minus_su;
    8485};
     
    105106  getColSolution = new double[model.getNumCols()];
    106107  CoinDisjointCopyN(model.getColSolution(), model.getNumCols(), getColSolution);
    107 
    108   getObjCoefficients = NULL;
     108  rowDualRay = NULL;
     109  colDualRay = NULL;
     110
    109111  yb_plus_rl_minus_su = 0;
    110112 
     
    121123  // that flag is correct, so we can test it.
    122124  if (isProvenPrimalInfeasible) {
    123     getObjCoefficients = new double[model.getNumCols()];
    124     CoinDisjointCopyN(model.getObjCoefficients(), model.getNumCols(), getObjCoefficients);
    125125    const std::vector<double*> dual_rays = model.getDualRays(1);
    126126    if (dual_rays.size() == 0) {
     
    128128      return;
    129129    }
    130     const double* dual_ray = dual_rays[0];
    131     const double direction = model.getObjSense();
    132130    const double* rub = model.getRowUpper();
    133131    const double* rlb = model.getRowLower();
    134132    const double* cub = model.getColUpper();
    135133    const double* clb = model.getColLower();
    136     const double* dj  = model.getReducedCost();
    137     const double* obj = model.getObjCoefficients();
    138134    const int numRows = model.getNumRows();
    139135    const int numCols = model.getNumCols();
     136    rowDualRay = dual_rays[0];
     137    colDualRay = new double[model.getNumCols()];
     138#if 1
     139    // WARNING  ==== WARNING ==== WARNING ==== WARNING ==== WARNING
     140    // For some reason clp returns the negative of the dual ray...
    140141    for (int i = 0; i < numRows; ++i) {
    141       const double ray_i = dual_ray[i];
    142       if (ray_i > 1e-6) {
     142      rowDualRay[i] = -rowDualRay[i];
     143    }
     144#endif
     145    const CoinPackedMatrix* m = model.getMatrixByCol();
     146    m->transposeTimes(rowDualRay, colDualRay);
     147   
     148    for (int i = 0; i < numRows; ++i) {
     149      const double ray_i = rowDualRay[i];
     150      if (ray_i > 1e-8) {
    143151        yb_plus_rl_minus_su += ray_i*rlb[i];
    144       } else if (ray_i < 1e-6) {
     152        assert(rlb[i] > -1e29);
     153      } else if (ray_i < 1e-8) {
    145154        yb_plus_rl_minus_su += ray_i*rub[i];
     155        assert(rub[i] < 1e29);
    146156      }
    147157    }
    148158    for (int j = 0; j < numCols; ++j) {
    149       const double yA_j = dj[j] - obj[j];
    150       if (direction * yA_j > 1e-6) {
    151         yb_plus_rl_minus_su -= yA_j*cub[j];
    152       } else if (direction * yA_j < -1e-6) {
    153         yb_plus_rl_minus_su -= yA_j*clb[j];
    154       }
    155     }
    156     for (int i = dual_rays.size()-1; i >= 0; --i) {
     159      if (colDualRay[j] > 1e-8) { //
     160        yb_plus_rl_minus_su -= colDualRay[j] * cub[j];
     161        assert(cub[j] < 1e29);
     162      } else if (colDualRay[j] < -1e-8) {
     163        yb_plus_rl_minus_su -= colDualRay[j] * clb[j];
     164        assert(clb[j] > -1e29);
     165      }
     166    }
     167    for (int i = dual_rays.size()-1; i > 0; --i) {
    157168      delete[] dual_rays[i];
    158169    }
     170    assert(yb_plus_rl_minus_su > 1e-8);
    159171  }
    160172}
     
    171183  delete[] getColUpper;
    172184  delete[] getColSolution;
    173   delete[] getObjCoefficients;
     185  delete[] colDualRay;
     186  delete[] rowDualRay;
    174187}
    175188
     
    226239  // Tests if we can switch this node (this is the parent) with its parent
    227240  bool canSwitchParentWithGrandparent(const int* which,
    228                                       const LPresult& lpres,
     241                                      LPresult& lpres,
    229242                                      const int * original_lower,
    230243                                      const int * original_upper,
     
    915928bool
    916929DBNodeSimple::canSwitchParentWithGrandparent(const int* which,
    917                                              const LPresult& lpres,
     930                                             LPresult& lpres,
    918931                                             const int * original_lower,
    919932                                             const int * original_upper,
     
    10251038  */
    10261039   
    1027   const double* dj = lpres.getReducedCost;
    1028   const double* obj = lpres.getObjCoefficients;
    1029   const double yA_x = dj[x] - obj[x];
    1030 
    1031   if (direction > 0) { // minimization problem
    1032     if (parent_is_down_child) {
    1033       const double s_x = CoinMax(yA_x, -1e-8);
    1034       if (s_x < 1e-6) {
    1035         return true;
    1036       }
    1037       if (lpres.yb_plus_rl_minus_su < 1e-8) {
    1038         printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
    1039         return false;
    1040       }
    1041       const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
    1042       const double u_x_without_GP = greatgrandparent_id >= 0 ?
    1043         branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
    1044         original_upper[GP_brvar];
    1045       return max_u_x > u_x_without_GP - 1e-8;
    1046     } else {
    1047       const double r_x = CoinMax(yA_x, -1e-8);
    1048       if (r_x < 1e-6) {
    1049         return true;
    1050       }
    1051       if (lpres.yb_plus_rl_minus_su < 1e-8) {
    1052         printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
    1053         return false;
    1054       }
    1055       const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
    1056       const double l_x_without_GP = greatgrandparent_id >= 0 ?
    1057         branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
    1058         original_lower[GP_brvar];
    1059       return min_l_x < l_x_without_GP + 1e-8;
    1060     }
    1061   } else { // maximization problem
    1062     if (parent_is_down_child) {
    1063       const double s_x = CoinMin(yA_x, 1e-8);
    1064       if (s_x > -1e-6) {
    1065         return true;
    1066       }
    1067       if (lpres.yb_plus_rl_minus_su > -1e-8) {
    1068         printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
    1069         return false;
    1070       }
    1071       const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
    1072       const double u_x_without_GP = greatgrandparent_id >= 0 ?
    1073         branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
    1074         original_upper[GP_brvar];
    1075       return max_u_x > u_x_without_GP - 1e-8;
    1076     } else {
    1077       const double r_x = CoinMin(yA_x, 1e-8);
    1078       if (r_x < -1e-6) {
    1079         return true;
    1080       }
    1081       if (lpres.yb_plus_rl_minus_su > -1e-8) {
    1082         printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
    1083         return false;
    1084       }
    1085       const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
    1086       const double l_x_without_GP = greatgrandparent_id >= 0 ?
    1087         branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
    1088         original_lower[GP_brvar];
    1089       return min_l_x < l_x_without_GP + 1e-8;
    1090     }
    1091   }
     1040  const double* colDualRay = lpres.colDualRay;
     1041  const double ub_without_GP = greatgrandparent_id >= 0 ?
     1042    branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
     1043    original_upper[GP_brvar];
     1044  const double lb_without_GP = greatgrandparent_id >= 0 ?
     1045    branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
     1046    original_lower[GP_brvar];
     1047
     1048  if (parent_is_down_child) { // upper bound on x is forced
     1049    if (colDualRay[x] < 1e-8) { // if yA_x] is non-positive then the var is
     1050                                // at its lower bound (or has 0 red cost)
     1051                                // so same dual ray still proves infeas
     1052      return true;
     1053    }
     1054    const double max_ub_increase = lpres.yb_plus_rl_minus_su / colDualRay[x];
     1055    assert(max_ub_increase >= 1e-8);
     1056    const bool willSwitch =
     1057      max_ub_increase > ub_without_GP - lpres.getColUpper[x] + 1e-8;
     1058    if (willSwitch) {
     1059      // The switch will happen so we might as well update the "infeasibility"
     1060      lpres.yb_plus_rl_minus_su -= colDualRay[x] * (ub_without_GP - lpres.getColUpper[x]);
     1061    }
     1062    return willSwitch;
     1063  } else { // lower bound on x is forced
     1064    if (colDualRay[x] > -1e-8) { // if yA_x is non-negative then the var is
     1065                                 // at its upper bound (or has 0 red cost)
     1066                                 // so same dual ray still proves infeas
     1067      return true;
     1068    }
     1069    const double max_lb_decrease = - lpres.yb_plus_rl_minus_su / colDualRay[x];
     1070    assert(max_lb_decrease >= 1e-8);
     1071    const bool willSwitch =
     1072      max_lb_decrease > lpres.getColLower[x] - lb_without_GP + 1e-8;
     1073    if (willSwitch) {
     1074      // The switch will happen so we might as well update the "infeasibility"
     1075      lpres.yb_plus_rl_minus_su += colDualRay[x] * (lpres.getColLower[x] - lb_without_GP);
     1076    }
     1077    return willSwitch;
     1078  }
     1079  // THINK: the above is definitely good for minimization problems. Is it good
     1080  // for max?
    10921081
    10931082  return true; // to placate some compilers
     
    16581647          // This is horribly looking... Get rid of it when properly
    16591648          // debugged...
     1649#if 1
     1650          const bool mustResolve =
     1651            model.isDualObjectiveLimitReached() && !model.isProvenOptimal();
     1652          double oldlimit = 0;
     1653
     1654          if (mustResolve) {
     1655            // THINK: Something faster would be better...
     1656            model.getDblParam(OsiDualObjectiveLimit, oldlimit);
     1657            model.setDblParam(OsiDualObjectiveLimit, 1e100);
     1658            model.resolve();
     1659          }
    16601660          LPresult lpres1(model);
     1661          if (mustResolve) {
     1662            lpres.isDualObjectiveLimitReached = true;
     1663            model.setDblParam(OsiDualObjectiveLimit, oldlimit);
     1664          }
    16611665          printf("After resolve: opt: %c   inf: %c   dual_bd: %c\n",
    16621666                 lpres1.isProvenOptimal ? 'T' : 'F',
     
    16721676                 (lpres.isProvenOptimal && model.isProvenOptimal() &&
    16731677                  lpres.getObjValue == model.getObjValue()));
     1678#endif
    16741679          printf("Finished moving node %d up by %i levels.\n", node.node_id_, cnt);
    16751680        }
Note: See TracChangeset for help on using the changeset viewer.