Changeset 1564
- Timestamp:
- Dec 30, 2010 12:45:15 PM (10 years ago)
- Location:
- trunk/Cbc/src
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Cbc/src/CbcHeuristic.hpp
r1505 r1564 249 249 inline void setWhereFrom(int value) { 250 250 whereFrom_ = value; 251 } 252 inline int whereFrom() const { 253 return whereFrom_; 251 254 } 252 255 /** Upto this depth we call the tree shallow and the heuristic can be called -
trunk/Cbc/src/CbcHeuristicGreedy.cpp
r1286 r1564 861 861 } 862 862 863 863 // Default Constructor 864 CbcHeuristicGreedySOS::CbcHeuristicGreedySOS() 865 : CbcHeuristic() 866 { 867 originalRhs_ = NULL; 868 // matrix will automatically be empty 869 originalNumberRows_ = 0; 870 algorithm_ = 0; 871 numberTimes_ = 100; 872 } 873 874 // Constructor from model 875 CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(CbcModel & model) 876 : CbcHeuristic(model) 877 { 878 gutsOfConstructor(&model); 879 algorithm_ = 2; 880 numberTimes_ = 100; 881 whereFrom_ = 1; 882 } 883 884 // Destructor 885 CbcHeuristicGreedySOS::~CbcHeuristicGreedySOS () 886 { 887 delete [] originalRhs_; 888 } 889 890 // Clone 891 CbcHeuristic * 892 CbcHeuristicGreedySOS::clone() const 893 { 894 return new CbcHeuristicGreedySOS(*this); 895 } 896 // Guts of constructor from a CbcModel 897 void 898 CbcHeuristicGreedySOS::gutsOfConstructor(CbcModel * model) 899 { 900 model_ = model; 901 // Get a copy of original matrix 902 assert(model->solver()); 903 if (model->solver()->getNumRows()) { 904 matrix_ = *model->solver()->getMatrixByCol(); 905 } 906 originalNumberRows_ = model->solver()->getNumRows(); 907 originalRhs_ = new double [originalNumberRows_]; 908 } 909 // Create C++ lines to get to current state 910 void 911 CbcHeuristicGreedySOS::generateCpp( FILE * fp) 912 { 913 CbcHeuristicGreedySOS other; 914 fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n"); 915 fprintf(fp, "3 CbcHeuristicGreedySOS heuristicGreedySOS(*cbcModel);\n"); 916 CbcHeuristic::generateCpp(fp, "heuristicGreedySOS"); 917 if (algorithm_ != other.algorithm_) 918 fprintf(fp, "3 heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_); 919 else 920 fprintf(fp, "4 heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_); 921 if (numberTimes_ != other.numberTimes_) 922 fprintf(fp, "3 heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_); 923 else 924 fprintf(fp, "4 heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_); 925 fprintf(fp, "3 cbcModel->addHeuristic(&heuristicGreedySOS);\n"); 926 } 927 928 // Copy constructor 929 CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(const CbcHeuristicGreedySOS & rhs) 930 : 931 CbcHeuristic(rhs), 932 matrix_(rhs.matrix_), 933 originalNumberRows_(rhs.originalNumberRows_), 934 algorithm_(rhs.algorithm_), 935 numberTimes_(rhs.numberTimes_) 936 { 937 originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_); 938 } 939 940 // Assignment operator 941 CbcHeuristicGreedySOS & 942 CbcHeuristicGreedySOS::operator=( const CbcHeuristicGreedySOS & rhs) 943 { 944 if (this != &rhs) { 945 CbcHeuristic::operator=(rhs); 946 matrix_ = rhs.matrix_; 947 originalNumberRows_ = rhs.originalNumberRows_; 948 algorithm_ = rhs.algorithm_; 949 numberTimes_ = rhs.numberTimes_; 950 delete [] originalRhs_; 951 originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_); 952 } 953 return *this; 954 } 955 // Returns 1 if solution, 0 if not 956 int 957 CbcHeuristicGreedySOS::solution(double & solutionValue, 958 double * betterSolution) 959 { 960 numCouldRun_++; 961 if (!model_) 962 return 0; 963 // See if to do 964 if (!when() || (when() == 1 && model_->phase() != 1)) 965 return 0; // switched off 966 if (model_->getNodeCount() > numberTimes_) 967 return 0; 968 // See if at root node 969 bool atRoot = model_->getNodeCount() == 0; 970 int passNumber = model_->getCurrentPassNumber(); 971 if (atRoot && passNumber != 1) 972 return 0; 973 OsiSolverInterface * solver = model_->solver(); 974 int numberColumns = solver->getNumCols(); 975 // This is number of rows when matrix was passed in 976 int numberRows = originalNumberRows_; 977 if (!numberRows) 978 return 0; // switched off 979 980 const double * columnLower = solver->getColLower(); 981 const double * columnUpper = solver->getColUpper(); 982 // modified rhs 983 double * rhs = CoinCopyOfArray(originalRhs_,numberRows); 984 // Column copy 985 const double * element = matrix_.getElements(); 986 const int * row = matrix_.getIndices(); 987 const CoinBigIndex * columnStart = matrix_.getVectorStarts(); 988 const int * columnLength = matrix_.getVectorLengths(); 989 // If bit set then use current 990 if ((algorithm_&1)!=0) { 991 const CoinPackedMatrix * matrix = solver->getMatrixByCol(); 992 element = matrix->getElements(); 993 row = matrix->getIndices(); 994 columnStart = matrix->getVectorStarts(); 995 columnLength = matrix->getVectorLengths(); 996 rhs = new double [numberRows]; 997 const double * rowLower = solver->getRowLower(); 998 const double * rowUpper = solver->getRowUpper(); 999 bool good = true; 1000 for (int iRow = 0; iRow < numberRows; iRow++) { 1001 if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) { 1002 // SOS 1003 rhs[iRow]=-1.0; 1004 } else if (rowLower[iRow] > 0.0) { 1005 good = false; 1006 } else if (rowUpper[iRow] < 0.0) { 1007 good = false; 1008 } else { 1009 rhs[iRow]=rowUpper[iRow]; 1010 } 1011 } 1012 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1013 if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0) 1014 good = false; 1015 CoinBigIndex j; 1016 int nSOS=0; 1017 if (!solver->isInteger(iColumn)) 1018 good = false; 1019 for (j = columnStart[iColumn]; 1020 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1021 if (element[j] < 0.0) 1022 good = false; 1023 int iRow = row[j]; 1024 if (rhs[iRow]==-1.0) { 1025 if (element[j] != 1.0) 1026 good = false; 1027 nSOS++; 1028 } 1029 } 1030 if (nSOS!=1) 1031 good = false; 1032 } 1033 if (!good) { 1034 delete [] rhs; 1035 setWhen(0); // switch off 1036 return 0; 1037 } 1038 } 1039 const double * solution = solver->getColSolution(); 1040 const double * objective = solver->getObjCoefficients(); 1041 double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance); 1042 double primalTolerance; 1043 solver->getDblParam(OsiPrimalTolerance, primalTolerance); 1044 1045 numRuns_++; 1046 assert (numberRows == matrix_.getNumRows()); 1047 int iRow, iColumn; 1048 double direction = solver->getObjSense(); 1049 double * slackCost = new double [numberRows]; 1050 double * modifiedCost = CoinCopyOfArray(objective,numberColumns); 1051 for (int iRow = 0;iRow < numberRows; iRow++) 1052 slackCost[iRow]=1.0e30; 1053 // Take off cost of gub slack 1054 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1055 if (columnLength[iColumn] == 1) { 1056 // SOS slack 1057 double cost = direction*objective[iColumn]; 1058 int iRow = row[columnStart[iColumn]]; 1059 assert (rhs[iRow]<0.0); 1060 slackCost[iRow]=CoinMin(slackCost[iRow],cost); 1061 } 1062 } 1063 double offset2 = 0.0; 1064 for (int iRow = 0;iRow < numberRows; iRow++) { 1065 if( slackCost[iRow] == 1.0e30) { 1066 slackCost[iRow]=0.0; 1067 } else { 1068 offset2 += slackCost[iRow]; 1069 rhs[iRow] = -2.0; 1070 } 1071 } 1072 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1073 double cost = direction * modifiedCost[iColumn]; 1074 CoinBigIndex j; 1075 for (j = columnStart[iColumn]; 1076 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1077 int iRow = row[j]; 1078 if (rhs[iRow]<0.0) { 1079 cost -= slackCost[iRow]; 1080 } 1081 } 1082 modifiedCost[iColumn] = cost; 1083 } 1084 delete [] slackCost; 1085 double offset; 1086 solver->getDblParam(OsiObjOffset, offset); 1087 double newSolutionValue = -offset+offset2; 1088 int returnCode = 0; 1089 1090 1091 // Get solution array for heuristic solution 1092 double * newSolution = new double [numberColumns]; 1093 double * rowActivity = new double[numberRows]; 1094 memset(rowActivity, 0, numberRows*sizeof(double)); 1095 if ((algorithm_&2)==0) { 1096 // get solution as small as possible 1097 for (iColumn = 0; iColumn < numberColumns; iColumn++) 1098 newSolution[iColumn] = columnLower[iColumn]; 1099 } else { 1100 // Get rounded down solution 1101 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1102 double value = solution[iColumn]; 1103 // Round down integer 1104 if (fabs(floor(value + 0.5) - value) < integerTolerance) { 1105 value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn])); 1106 } else { 1107 value = CoinMax(floor(value), columnLower[iColumn]); 1108 } 1109 // make sure clean 1110 value = CoinMin(value, columnUpper[iColumn]); 1111 value = CoinMax(value, columnLower[iColumn]); 1112 newSolution[iColumn] = value; 1113 } 1114 } 1115 double * contribution = new double [numberColumns]; 1116 int * which = new int [numberColumns]; 1117 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1118 CoinBigIndex j; 1119 double value = newSolution[iColumn]; 1120 double cost = modifiedCost[iColumn]; 1121 double forSort = 0.0; 1122 bool hasSlack=false; 1123 newSolutionValue += value * cost; 1124 for (j = columnStart[iColumn]; 1125 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1126 int iRow = row[j]; 1127 rowActivity[iRow] += value * element[j]; 1128 if (rhs[iRow] >0.0) { 1129 forSort += element[j]; 1130 } else if (rhs[iRow] == -2.0) { 1131 hasSlack = true; 1132 } 1133 } 1134 if (forSort && value == 0.0 && columnUpper[iColumn]) { 1135 if (hasSlack) { 1136 if (cost>=0.0) { 1137 forSort = 2.0e30; 1138 } else { 1139 forSort = cost/forSort; 1140 } 1141 } else { 1142 forSort = cost/forSort; 1143 } 1144 } else { 1145 // put at end 1146 forSort = 1.0e30; 1147 } 1148 which[iColumn]=iColumn; 1149 contribution[iColumn]= forSort; 1150 } 1151 CoinSort_2(contribution,contribution+numberColumns,which); 1152 // Go through columns 1153 for (int jColumn = 0; jColumn < numberColumns; jColumn++) { 1154 int iColumn = which[jColumn]; 1155 double value = newSolution[iColumn]; 1156 if (value) 1157 continue; 1158 bool possible = true; 1159 CoinBigIndex j; 1160 for (j = columnStart[iColumn]; 1161 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1162 int iRow = row[j]; 1163 if (rhs[iRow]<0.0&&rowActivity[iRow]) { 1164 possible = false; 1165 } else if (rhs[iRow]>=0.0) { 1166 double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8; 1167 if (gap<element[j]) 1168 possible = false; 1169 } 1170 } 1171 if (possible) { 1172 // Increase chosen column 1173 newSolution[iColumn] = 1.0; 1174 double cost = modifiedCost[iColumn]; 1175 newSolutionValue += cost; 1176 for (CoinBigIndex j = columnStart[iColumn]; 1177 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1178 int iRow = row[j]; 1179 rowActivity[iRow] += element[j]; 1180 } 1181 } 1182 } 1183 if (newSolutionValue < solutionValue) { 1184 // check feasible 1185 const double * rowLower = solver->getRowLower(); 1186 const double * rowUpper = solver->getRowUpper(); 1187 memset(rowActivity, 0, numberRows*sizeof(double)); 1188 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1189 CoinBigIndex j; 1190 double value = newSolution[iColumn]; 1191 if (value) { 1192 for (j = columnStart[iColumn]; 1193 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1194 int iRow = row[j]; 1195 rowActivity[iRow] += value * element[j]; 1196 } 1197 } 1198 } 1199 // check was approximately feasible 1200 bool feasible = true; 1201 for (iRow = 0; iRow < numberRows; iRow++) { 1202 if (rowActivity[iRow] < rowLower[iRow]) { 1203 if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) 1204 feasible = false; 1205 } else if (rowActivity[iRow] > rowUpper[iRow]) { 1206 if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) 1207 feasible = false; 1208 } 1209 } 1210 if (feasible) { 1211 // new solution 1212 memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); 1213 solutionValue = newSolutionValue; 1214 //printf("** Solution of %g found by rounding\n",newSolutionValue); 1215 returnCode = 1; 1216 } else { 1217 // Can easily happen 1218 //printf("Debug CbcHeuristicGreedySOS giving bad solution\n"); 1219 } 1220 } 1221 delete [] newSolution; 1222 delete [] rowActivity; 1223 delete [] modifiedCost; 1224 delete [] contribution; 1225 delete [] which; 1226 delete [] rhs; 1227 return returnCode; 1228 } 1229 // update model 1230 void CbcHeuristicGreedySOS::setModel(CbcModel * model) 1231 { 1232 delete [] originalRhs_; 1233 gutsOfConstructor(model); 1234 validate(); 1235 } 1236 // Resets stuff if model changes 1237 void 1238 CbcHeuristicGreedySOS::resetModel(CbcModel * model) 1239 { 1240 delete [] originalRhs_; 1241 gutsOfConstructor(model); 1242 } 1243 // Validate model i.e. sets when_ to 0 if necessary (may be NULL) 1244 void 1245 CbcHeuristicGreedySOS::validate() 1246 { 1247 if (model_ && when() < 10) { 1248 if (model_->numberIntegers() != 1249 model_->numberObjects() && (model_->numberObjects() || 1250 (model_->specialOptions()&1024) == 0)) { 1251 int numberOdd = 0; 1252 for (int i = 0; i < model_->numberObjects(); i++) { 1253 if (!model_->object(i)->canDoHeuristics()) 1254 numberOdd++; 1255 } 1256 if (numberOdd) 1257 setWhen(0); 1258 } 1259 // Only works if coefficients positive and all rows L or SOS 1260 OsiSolverInterface * solver = model_->solver(); 1261 const double * columnUpper = solver->getColUpper(); 1262 const double * columnLower = solver->getColLower(); 1263 const double * rowLower = solver->getRowLower(); 1264 const double * rowUpper = solver->getRowUpper(); 1265 1266 int numberRows = solver->getNumRows(); 1267 // Column copy 1268 const double * element = matrix_.getElements(); 1269 const int * row = matrix_.getIndices(); 1270 const CoinBigIndex * columnStart = matrix_.getVectorStarts(); 1271 const int * columnLength = matrix_.getVectorLengths(); 1272 bool good = true; 1273 assert (originalRhs_); 1274 for (int iRow = 0; iRow < numberRows; iRow++) { 1275 if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) { 1276 // SOS 1277 originalRhs_[iRow]=-1.0; 1278 } else if (rowLower[iRow] > 0.0) { 1279 good = false; 1280 } else if (rowUpper[iRow] < 0.0) { 1281 good = false; 1282 } else { 1283 originalRhs_[iRow]=rowUpper[iRow]; 1284 } 1285 } 1286 int numberColumns = solver->getNumCols(); 1287 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1288 if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0) 1289 good = false; 1290 CoinBigIndex j; 1291 int nSOS=0; 1292 if (!solver->isInteger(iColumn)) 1293 good = false; 1294 for (j = columnStart[iColumn]; 1295 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1296 if (element[j] < 0.0) 1297 good = false; 1298 int iRow = row[j]; 1299 if (originalRhs_[iRow]==-1.0) { 1300 if (element[j] != 1.0) 1301 good = false; 1302 nSOS++; 1303 } 1304 } 1305 if (nSOS!=1) 1306 good = false; 1307 } 1308 if (!good) 1309 setWhen(0); // switch off 1310 } 1311 } -
trunk/Cbc/src/CbcHeuristicGreedy.hpp
r1432 r1564 186 186 }; 187 187 188 /** Greedy heuristic for SOS and L rows (and positive elements) 189 */ 190 191 class CbcHeuristicGreedySOS : public CbcHeuristic { 192 public: 193 194 // Default Constructor 195 CbcHeuristicGreedySOS (); 196 197 /* Constructor with model - assumed before cuts 198 Initial version does not do Lps 199 */ 200 CbcHeuristicGreedySOS (CbcModel & model); 201 202 // Copy constructor 203 CbcHeuristicGreedySOS ( const CbcHeuristicGreedySOS &); 204 205 // Destructor 206 ~CbcHeuristicGreedySOS (); 207 208 /// Clone 209 virtual CbcHeuristic * clone() const; 210 /// Assignment operator 211 CbcHeuristicGreedySOS & operator=(const CbcHeuristicGreedySOS& rhs); 212 /// Create C++ lines to get to current state 213 virtual void generateCpp( FILE * fp) ; 214 215 /// update model (This is needed if cliques update matrix etc) 216 virtual void setModel(CbcModel * model); 217 218 using CbcHeuristic::solution ; 219 /** returns 0 if no solution, 1 if valid solution. 220 Sets solution values if good, sets objective value (only if good) 221 We leave all variables which are at one at this node of the 222 tree to that value and will 223 initially set all others to zero. We then sort all variables in order of their cost 224 divided by the number of entries in rows which are not yet covered. We randomize that 225 value a bit so that ties will be broken in different ways on different runs of the heuristic. 226 We then choose the best one and set it to one and repeat the exercise. 227 228 */ 229 virtual int solution(double & objectiveValue, 230 double * newSolution); 231 /// Validate model i.e. sets when_ to 0 if necessary (may be NULL) 232 virtual void validate() ; 233 /// Resets stuff if model changes 234 virtual void resetModel(CbcModel * model); 235 /* Algorithm 236 Bits 237 1 bit - use current model, otherwise original 238 2 - use current solution as starting point, otherwise pure greedy 239 4 - use duals to modify greedy 240 8 - use duals on GUB/SOS in special way 241 */ 242 inline int algorithm() const { 243 return algorithm_; 244 } 245 inline void setAlgorithm(int value) { 246 algorithm_ = value; 247 } 248 // Only do this many times 249 inline int numberTimes() const { 250 return numberTimes_; 251 } 252 inline void setNumberTimes(int value) { 253 numberTimes_ = value; 254 } 255 256 protected: 257 /// Guts of constructor from a CbcModel 258 void gutsOfConstructor(CbcModel * model); 259 // Data 260 261 // Original RHS - if -1.0 then SOS otherwise <= value 262 double * originalRhs_; 263 // Original matrix by column 264 CoinPackedMatrix matrix_; 265 // original number of rows 266 int originalNumberRows_; 267 /* Algorithm 268 */ 269 int algorithm_; 270 /// Do this many times 271 int numberTimes_; 272 273 }; 274 188 275 189 276 #endif -
trunk/Cbc/src/CbcHeuristicRENS.cpp
r1499 r1564 14 14 #include "CbcHeuristicRENS.hpp" 15 15 #include "CoinWarmStartBasis.hpp" 16 #include "CoinSort.hpp" 16 17 #include "CbcBranchActual.hpp" 17 18 #include "CbcStrategy.hpp" … … 89 90 90 91 const double * currentSolution = solver->getColSolution(); 91 const double * dj = solver->getReducedCost(); 92 int numberColumns = solver->getNumCols(); 93 double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns); 92 94 OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone(); 93 int numberColumns = newSolver->getNumCols();94 95 double direction = newSolver->getObjSense(); 95 96 int type = rensType_&15; … … 97 98 const double * colLower = newSolver->getColLower(); 98 99 const double * colUpper = newSolver->getColUpper(); 99 if ((type&3)==3) { 100 double total=0.0; 101 int n=0; 102 CoinWarmStartBasis * basis = 103 dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; 104 if (basis) { 105 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 106 if (colUpper[iColumn]>colLower[iColumn]&& 107 basis->getStructStatus(iColumn) != 108 CoinWarmStartBasis::basic) { 109 n++; 110 total += fabs(dj[iColumn]); 100 int numberFixed = 0; 101 if ((type&7)==3) { 102 double total=0.0; 103 int n=0; 104 CoinWarmStartBasis * basis = 105 dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; 106 if (basis&&basis->getNumArtificial()) { 107 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 108 if (colUpper[iColumn]>colLower[iColumn]&& 109 basis->getStructStatus(iColumn) != 110 CoinWarmStartBasis::basic) { 111 n++; 112 total += fabs(dj[iColumn]); 113 } 114 } 115 if (n) 116 djTolerance = (0.01*total)/static_cast<double>(n); 117 delete basis; 118 } 119 } else if ((type&7)==4) { 120 double * sort = new double [numberColumns]; 121 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 122 sort[iColumn]=1.0e30; 123 if (colUpper[iColumn]>colLower[iColumn]) { 124 sort[iColumn] = fabs(dj[iColumn]); 125 } 126 } 127 std::sort(sort,sort+numberColumns); 128 int last = static_cast<int>(numberColumns*fractionSmall_); 129 djTolerance = sort[last]; 130 delete [] sort; 131 } else if ((type&7)==5||(type&7)==6) { 132 // SOS type fixing 133 bool fixSets = (type&7)==5; 134 CoinWarmStartBasis * basis = 135 dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; 136 if (basis&&basis->getNumArtificial()) { 137 //const double * rowLower = solver->getRowLower(); 138 const double * rowUpper = solver->getRowUpper(); 139 140 int numberRows = solver->getNumRows(); 141 // Column copy 142 const CoinPackedMatrix * matrix = solver->getMatrixByCol(); 143 const double * element = matrix->getElements(); 144 const int * row = matrix->getIndices(); 145 const CoinBigIndex * columnStart = matrix->getVectorStarts(); 146 const int * columnLength = matrix->getVectorLengths(); 147 double * bestDj = new double [numberRows]; 148 for (int i=0;i<numberRows;i++) { 149 if (rowUpper[i]==1.0) 150 bestDj[i]=1.0e20; 151 else 152 bestDj[i]=1.0e30; 153 } 154 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 155 if (colUpper[iColumn]>colLower[iColumn]) { 156 CoinBigIndex j; 157 if (currentSolution[iColumn]>1.0e-6&& 158 currentSolution[iColumn]<0.999999) { 159 for (j = columnStart[iColumn]; 160 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 161 int iRow = row[j]; 162 bestDj[iRow]=1.0e30; 163 } 164 } else if ( basis->getStructStatus(iColumn) != 165 CoinWarmStartBasis::basic) { 166 for (j = columnStart[iColumn]; 167 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 168 int iRow = row[j]; 169 if (bestDj[iRow]<1.0e30) { 170 if (element[j] != 1.0) 171 bestDj[iRow]=1.0e30; 172 else 173 bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]); 174 } 111 175 } 112 176 } 113 if (n) 114 djTolerance = (0.01*total)/static_cast<double>(n); 115 delete basis; 116 } 177 } 178 } 179 int nSOS=0; 180 double * sort = new double [numberRows]; 181 for (int i=0;i<numberRows;i++) { 182 if (bestDj[i]<1.0e30) { 183 sort[nSOS++]=bestDj[i]; 184 } 185 } 186 if (10*nSOS>8*numberRows) { 187 std::sort(sort,sort+nSOS); 188 int last = static_cast<int>(nSOS*0.9*fractionSmall_); 189 double tolerance = sort[last]; 190 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 191 if (colUpper[iColumn]>colLower[iColumn]) { 192 CoinBigIndex j; 193 if (currentSolution[iColumn]<=1.0e-6|| 194 currentSolution[iColumn]>=0.999999) { 195 if (fixSets) { 196 for (j = columnStart[iColumn]; 197 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 198 int iRow = row[j]; 199 if (bestDj[iRow]<1.0e30&&bestDj[iRow]>=tolerance) { 200 numberFixed++; 201 if (currentSolution[iColumn]<=1.0e-6) 202 newSolver->setColUpper(iColumn,0.0); 203 else if (currentSolution[iColumn]>=0.999999) 204 newSolver->setColLower(iColumn,1.0); 205 } 206 } 207 } else if (columnLength[iColumn]==1) { 208 // leave more slacks 209 int iRow = row[columnStart[iColumn]]; 210 if (bestDj[iRow]<1.0e30) { 211 // fake dj 212 dj[iColumn] *= 0.000001; 213 } 214 } 215 } 216 } 217 } 218 if (fixSets) 219 djTolerance = 1.0e30; 220 else 221 djTolerance = tolerance; 222 } 223 delete basis; 224 delete [] sort; 225 delete [] bestDj; 226 } 117 227 } 118 228 119 229 double primalTolerance; 120 230 solver->getDblParam(OsiPrimalTolerance, primalTolerance); 121 231 122 232 int i; 123 int numberFixed = 0;124 233 int numberTightened = 0; 125 234 int numberAtBound = 0; … … 183 292 #endif 184 293 } 294 delete [] dj; 185 295 if (numberFixed > numberIntegers / 5) { 186 296 if ( numberFixed < numberColumns / 5) {
Note: See TracChangeset
for help on using the changeset viewer.