 Timestamp:
 Jun 29, 2008 6:49:32 AM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/dynamicbranching/dynamicbranching.cpp
r994 r995 79 79 double* getColLower; 80 80 double* getColUpper; 81 double* rowDualRay; 82 double* colDualRay; 81 83 double* getColSolution; // FIXME Not needed, just for debugging 82 double* getObjCoefficients;83 84 double yb_plus_rl_minus_su; 84 85 }; … … 105 106 getColSolution = new double[model.getNumCols()]; 106 107 CoinDisjointCopyN(model.getColSolution(), model.getNumCols(), getColSolution); 107 108 getObjCoefficients = NULL; 108 rowDualRay = NULL; 109 colDualRay = NULL; 110 109 111 yb_plus_rl_minus_su = 0; 110 112 … … 121 123 // that flag is correct, so we can test it. 122 124 if (isProvenPrimalInfeasible) { 123 getObjCoefficients = new double[model.getNumCols()];124 CoinDisjointCopyN(model.getObjCoefficients(), model.getNumCols(), getObjCoefficients);125 125 const std::vector<double*> dual_rays = model.getDualRays(1); 126 126 if (dual_rays.size() == 0) { … … 128 128 return; 129 129 } 130 const double* dual_ray = dual_rays[0];131 const double direction = model.getObjSense();132 130 const double* rub = model.getRowUpper(); 133 131 const double* rlb = model.getRowLower(); 134 132 const double* cub = model.getColUpper(); 135 133 const double* clb = model.getColLower(); 136 const double* dj = model.getReducedCost();137 const double* obj = model.getObjCoefficients();138 134 const int numRows = model.getNumRows(); 139 135 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... 140 141 for (int i = 0; i < numRows; ++i) { 141 const double ray_i = dual_ray[i]; 142 if (ray_i > 1e6) { 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 > 1e8) { 143 151 yb_plus_rl_minus_su += ray_i*rlb[i]; 144 } else if (ray_i < 1e6) { 152 assert(rlb[i] > 1e29); 153 } else if (ray_i < 1e8) { 145 154 yb_plus_rl_minus_su += ray_i*rub[i]; 155 assert(rub[i] < 1e29); 146 156 } 147 157 } 148 158 for (int j = 0; j < numCols; ++j) { 149 const double yA_j = dj[j]  obj[j]; 150 if (direction * yA_j > 1e6) { 151 yb_plus_rl_minus_su = yA_j*cub[j]; 152 } else if (direction * yA_j < 1e6) { 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] > 1e8) { // 160 yb_plus_rl_minus_su = colDualRay[j] * cub[j]; 161 assert(cub[j] < 1e29); 162 } else if (colDualRay[j] < 1e8) { 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) { 157 168 delete[] dual_rays[i]; 158 169 } 170 assert(yb_plus_rl_minus_su > 1e8); 159 171 } 160 172 } … … 171 183 delete[] getColUpper; 172 184 delete[] getColSolution; 173 delete[] getObjCoefficients; 185 delete[] colDualRay; 186 delete[] rowDualRay; 174 187 } 175 188 … … 226 239 // Tests if we can switch this node (this is the parent) with its parent 227 240 bool canSwitchParentWithGrandparent(const int* which, 228 constLPresult& lpres,241 LPresult& lpres, 229 242 const int * original_lower, 230 243 const int * original_upper, … … 915 928 bool 916 929 DBNodeSimple::canSwitchParentWithGrandparent(const int* which, 917 constLPresult& lpres,930 LPresult& lpres, 918 931 const int * original_lower, 919 932 const int * original_upper, … … 1025 1038 */ 1026 1039 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, 1e8); 1034 if (s_x < 1e6) { 1035 return true; 1036 } 1037 if (lpres.yb_plus_rl_minus_su < 1e8) { 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  1e8; 1046 } else { 1047 const double r_x = CoinMax(yA_x, 1e8); 1048 if (r_x < 1e6) { 1049 return true; 1050 } 1051 if (lpres.yb_plus_rl_minus_su < 1e8) { 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 + 1e8; 1060 } 1061 } else { // maximization problem 1062 if (parent_is_down_child) { 1063 const double s_x = CoinMin(yA_x, 1e8); 1064 if (s_x > 1e6) { 1065 return true; 1066 } 1067 if (lpres.yb_plus_rl_minus_su > 1e8) { 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  1e8; 1076 } else { 1077 const double r_x = CoinMin(yA_x, 1e8); 1078 if (r_x < 1e6) { 1079 return true; 1080 } 1081 if (lpres.yb_plus_rl_minus_su > 1e8) { 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 + 1e8; 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] < 1e8) { // if yA_x] is nonpositive 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 >= 1e8); 1056 const bool willSwitch = 1057 max_ub_increase > ub_without_GP  lpres.getColUpper[x] + 1e8; 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] > 1e8) { // if yA_x is nonnegative 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 >= 1e8); 1071 const bool willSwitch = 1072 max_lb_decrease > lpres.getColLower[x]  lb_without_GP + 1e8; 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? 1092 1081 1093 1082 return true; // to placate some compilers … … 1658 1647 // This is horribly looking... Get rid of it when properly 1659 1648 // 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 } 1660 1660 LPresult lpres1(model); 1661 if (mustResolve) { 1662 lpres.isDualObjectiveLimitReached = true; 1663 model.setDblParam(OsiDualObjectiveLimit, oldlimit); 1664 } 1661 1665 printf("After resolve: opt: %c inf: %c dual_bd: %c\n", 1662 1666 lpres1.isProvenOptimal ? 'T' : 'F', … … 1672 1676 (lpres.isProvenOptimal && model.isProvenOptimal() && 1673 1677 lpres.getObjValue == model.getObjValue())); 1678 #endif 1674 1679 printf("Finished moving node %d up by %i levels.\n", node.node_id_, cnt); 1675 1680 }
Note: See TracChangeset
for help on using the changeset viewer.