Changeset 1564 for trunk/Cbc/src/CbcHeuristicGreedy.cpp
 Timestamp:
 Dec 30, 2010 12:45:15 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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.0e3, 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.0e8; 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 }
Note: See TracChangeset
for help on using the changeset viewer.