1 | /* $Id: Idiot.cpp 1753 2011-06-19 16:27:26Z stefan $ */ |
2 | // Copyright (C) 2002, International Business Machines |
3 | // Corporation and others. All Rights Reserved. |
4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
5 | |
6 | #include "CoinPragma.hpp" |
7 | #include <stdio.h> |
8 | #include <stdarg.h> |
9 | #include <stdlib.h> |
10 | #include <math.h> |
11 | #include "ClpPresolve.hpp" |
12 | #include "Idiot.hpp" |
13 | #include "CoinTime.hpp" |
14 | #include "CoinSort.hpp" |
15 | #include "CoinMessageHandler.hpp" |
16 | #include "CoinHelperFunctions.hpp" |
17 | // Redefine stuff for Clp |
18 | #ifndef OSI_IDIOT |
19 | #include "ClpMessage.hpp" |
20 | #define OsiObjOffset ClpObjOffset |
21 | #endif |
22 | /**** strategy 4 - drop, exitDrop and djTolerance all relative: |
23 | For first two major iterations these are small. Then: |
24 | |
25 | drop - exit a major iteration if drop over 5*checkFrequency < this is |
26 | used as info->drop*(10.0+fabs(last weighted objective)) |
27 | |
28 | exitDrop - exit idiot if feasible and drop < this is |
29 | used as info->exitDrop*(10.0+fabs(last objective)) |
30 | |
31 | djExit - exit a major iteration if largest dj (averaged over 5 checks) |
32 | drops below this - used as info->djTolerance*(10.0+fabs(last weighted objective) |
33 | |
34 | djFlag - mostly skip variables with bad dj worse than this => 2*djExit |
35 | |
36 | djTol - only look at variables with dj better than this => 0.01*djExit |
37 | ****************/ |
38 | |
39 | #define IDIOT_FIX_TOLERANCE 1e-6 |
40 | #define SMALL_IDIOT_FIX_TOLERANCE 1e-10 |
41 | int |
42 | Idiot::dropping(IdiotResult result, |
43 | double tolerance, |
44 | double small, |
45 | int *nbad) |
46 | { |
47 | if (result.infeas <= small) { |
48 | double value = CoinMax(fabs(result.objval), fabs(result.dropThis)) + 1.0; |
49 | if (result.dropThis > tolerance * value) { |
50 | *nbad = 0; |
51 | return 1; |
52 | } else { |
53 | (*nbad)++; |
54 | if (*nbad > 4) { |
55 | return 0; |
56 | } else { |
57 | return 1; |
58 | } |
59 | } |
60 | } else { |
61 | *nbad = 0; |
62 | return 1; |
63 | } |
64 | } |
65 | // Deals with whenUsed and slacks |
66 | int |
67 | Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd, |
68 | double * colsol, const double * lower, const double * upper, |
69 | const double * rowLower, const double * rowUpper, |
70 | const double * cost, const double * element, double fixTolerance, |
71 | double & objValue, double & infValue) |
72 | { |
73 | int n = 0; |
74 | if ((strategy_ & 16384) == 0) { |
75 | for (int i = ordinaryStart; i < ordinaryEnd; i++) { |
76 | if (colsol[i] > lower[i] + fixTolerance) { |
77 | if (colsol[i] < upper[i] - fixTolerance) { |
78 | n++; |
79 | } else { |
80 | colsol[i] = upper[i]; |
81 | } |
82 | whenUsed_[i] = iteration; |
83 | } else { |
84 | colsol[i] = lower[i]; |
85 | } |
86 | } |
87 | return n; |
88 | } else { |
89 | #ifdef COIN_DEVELOP |
90 | printf("entering inf %g, obj %g\n" , infValue, objValue); |
91 | #endif |
92 | int nrows = model_->getNumRows(); |
93 | int ncols = model_->getNumCols(); |
94 | int * posSlack = whenUsed_ + ncols; |
95 | int * negSlack = posSlack + nrows; |
96 | int * nextSlack = negSlack + nrows; |
97 | double * rowsol = reinterpret_cast<double *> (nextSlack + ncols); |
98 | memset(rowsol, 0, nrows * sizeof(double)); |
99 | #ifdef OSI_IDIOT |
100 | const CoinPackedMatrix * matrix = model_->getMatrixByCol(); |
101 | #else |
102 | ClpMatrixBase * matrix = model_->clpMatrix(); |
103 | #endif |
104 | const int * row = matrix->getIndices(); |
105 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
106 | const int * columnLength = matrix->getVectorLengths(); |
107 | //const double * element = matrix->getElements(); |
108 | int i; |
109 | objValue = 0.0; |
110 | infValue = 0.0; |
111 | for ( i = 0; i < ncols; i++) { |
112 | if (nextSlack[i] == -1) { |
113 | // not a slack |
114 | if (colsol[i] > lower[i] + fixTolerance) { |
115 | if (colsol[i] < upper[i] - fixTolerance) { |
116 | n++; |
117 | whenUsed_[i] = iteration; |
118 | } else { |
119 | colsol[i] = upper[i]; |
120 | } |
121 | whenUsed_[i] = iteration; |
122 | } else { |
123 | colsol[i] = lower[i]; |
124 | } |
125 | double value = colsol[i]; |
126 | if (value) { |
127 | objValue += cost[i] * value; |
128 | CoinBigIndex j; |
129 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
130 | int iRow = row[j]; |
131 | rowsol[iRow] += value * element[j]; |
132 | } |
133 | } |
134 | } |
135 | } |
136 | // temp fix for infinite lbs - just limit to -1000 |
137 | for (i = 0; i < nrows; i++) { |
138 | double rowSave = rowsol[i]; |
139 | int iCol; |
140 | iCol = posSlack[i]; |
141 | if (iCol >= 0) { |
142 | // slide all slack down |
143 | double rowValue = rowsol[i]; |
144 | CoinBigIndex j = columnStart[iCol]; |
145 | double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]); |
146 | rowSave += (colsol[iCol] - lowerValue) * element[j]; |
147 | colsol[iCol] = lowerValue; |
148 | while (nextSlack[iCol] >= 0) { |
149 | iCol = nextSlack[iCol]; |
150 | double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]); |
151 | j = columnStart[iCol]; |
152 | rowSave += (colsol[iCol] - lowerValue) * element[j]; |
153 | colsol[iCol] = lowerValue; |
154 | } |
155 | iCol = posSlack[i]; |
156 | while (rowValue < rowLower[i] && iCol >= 0) { |
157 | // want to increase |
158 | double distance = rowLower[i] - rowValue; |
159 | double value = element[columnStart[iCol]]; |
160 | double thisCost = cost[iCol]; |
161 | if (distance <= value*(upper[iCol] - colsol[iCol])) { |
162 | // can get there |
163 | double movement = distance / value; |
164 | objValue += movement * thisCost; |
165 | rowValue = rowLower[i]; |
166 | colsol[iCol] += movement; |
167 | } else { |
168 | // can't get there |
169 | double movement = upper[iCol] - colsol[iCol]; |
170 | objValue += movement * thisCost; |
171 | rowValue += movement * value; |
172 | colsol[iCol] = upper[iCol]; |
173 | iCol = nextSlack[iCol]; |
174 | } |
175 | } |
176 | if (iCol >= 0) { |
177 | // may want to carry on - because of cost? |
178 | while (iCol >= 0 && cost[iCol] < 0 && rowValue < rowUpper[i]) { |
179 | // want to increase |
180 | double distance = rowUpper[i] - rowValue; |
181 | double value = element[columnStart[iCol]]; |
182 | double thisCost = cost[iCol]; |
183 | if (distance <= value*(upper[iCol] - colsol[iCol])) { |
184 | // can get there |
185 | double movement = distance / value; |
186 | objValue += movement * thisCost; |
187 | rowValue = rowUpper[i]; |
188 | colsol[iCol] += movement; |
189 | iCol = -1; |
190 | } else { |
191 | // can't get there |
192 | double movement = upper[iCol] - colsol[iCol]; |
193 | objValue += movement * thisCost; |
194 | rowValue += movement * value; |
195 | colsol[iCol] = upper[iCol]; |
196 | iCol = nextSlack[iCol]; |
197 | } |
198 | } |
199 | if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance && |
200 | colsol[iCol] < upper[iCol] - fixTolerance) { |
201 | whenUsed_[i] = iteration; |
202 | n++; |
203 | } |
204 | } |
205 | rowsol[i] = rowValue; |
206 | } |
207 | iCol = negSlack[i]; |
208 | if (iCol >= 0) { |
209 | // slide all slack down |
210 | double rowValue = rowsol[i]; |
211 | CoinBigIndex j = columnStart[iCol]; |
212 | rowSave += (colsol[iCol] - lower[iCol]) * element[j]; |
213 | colsol[iCol] = lower[iCol]; |
214 | assert (lower[iCol] > -1.0e20); |
215 | while (nextSlack[iCol] >= 0) { |
216 | iCol = nextSlack[iCol]; |
217 | j = columnStart[iCol]; |
218 | rowSave += (colsol[iCol] - lower[iCol]) * element[j]; |
219 | colsol[iCol] = lower[iCol]; |
220 | } |
221 | iCol = negSlack[i]; |
222 | while (rowValue > rowUpper[i] && iCol >= 0) { |
223 | // want to increase |
224 | double distance = -(rowUpper[i] - rowValue); |
225 | double value = -element[columnStart[iCol]]; |
226 | double thisCost = cost[iCol]; |
227 | if (distance <= value*(upper[iCol] - lower[iCol])) { |
228 | // can get there |
229 | double movement = distance / value; |
230 | objValue += movement * thisCost; |
231 | rowValue = rowUpper[i]; |
232 | colsol[iCol] += movement; |
233 | } else { |
234 | // can't get there |
235 | double movement = upper[iCol] - lower[iCol]; |
236 | objValue += movement * thisCost; |
237 | rowValue -= movement * value; |
238 | colsol[iCol] = upper[iCol]; |
239 | iCol = nextSlack[iCol]; |
240 | } |
241 | } |
242 | if (iCol >= 0) { |
243 | // may want to carry on - because of cost? |
244 | while (iCol >= 0 && cost[iCol] < 0 && rowValue > rowLower[i]) { |
245 | // want to increase |
246 | double distance = -(rowLower[i] - rowValue); |
247 | double value = -element[columnStart[iCol]]; |
248 | double thisCost = cost[iCol]; |
249 | if (distance <= value*(upper[iCol] - colsol[iCol])) { |
250 | // can get there |
251 | double movement = distance / value; |
252 | objValue += movement * thisCost; |
253 | rowValue = rowLower[i]; |
254 | colsol[iCol] += movement; |
255 | iCol = -1; |
256 | } else { |
257 | // can't get there |
258 | double movement = upper[iCol] - colsol[iCol]; |
259 | objValue += movement * thisCost; |
260 | rowValue -= movement * value; |
261 | colsol[iCol] = upper[iCol]; |
262 | iCol = nextSlack[iCol]; |
263 | } |
264 | } |
265 | if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance && |
266 | colsol[iCol] < upper[iCol] - fixTolerance) { |
267 | whenUsed_[i] = iteration; |
268 | n++; |
269 | } |
270 | } |
271 | rowsol[i] = rowValue; |
272 | } |
273 | infValue += CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]); |
274 | // just change |
275 | rowsol[i] -= rowSave; |
276 | } |
277 | return n; |
278 | } |
279 | } |
280 | |
281 | /* returns -1 if none or start of costed slacks or -2 if |
282 | there are costed slacks but it is messy */ |
283 | static int countCostedSlacks(OsiSolverInterface * model) |
284 | { |
285 | #ifdef OSI_IDIOT |
286 | const CoinPackedMatrix * matrix = model->getMatrixByCol(); |
287 | #else |
288 | ClpMatrixBase * matrix = model->clpMatrix(); |
289 | #endif |
290 | const int * row = matrix->getIndices(); |
291 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
292 | const int * columnLength = matrix->getVectorLengths(); |
293 | const double * element = matrix->getElements(); |
294 | const double * rowupper = model->getRowUpper(); |
295 | int nrows = model->getNumRows(); |
296 | int ncols = model->getNumCols(); |
297 | int slackStart = ncols - nrows; |
298 | int nSlacks = nrows; |
299 | int i; |
300 | |
301 | if (ncols <= nrows) return -1; |
302 | while (1) { |
303 | for (i = 0; i < nrows; i++) { |
304 | int j = i + slackStart; |
305 | CoinBigIndex k = columnStart[j]; |
306 | if (columnLength[j] == 1) { |
307 | if (row[k] != i || element[k] != 1.0) { |
308 | nSlacks = 0; |
309 | break; |
310 | } |
311 | } else { |
312 | nSlacks = 0; |
313 | break; |
314 | } |
315 | if (rowupper[i] <= 0.0) { |
316 | nSlacks = 0; |
317 | break; |
318 | } |
319 | } |
320 | if (nSlacks || !slackStart) break; |
321 | slackStart = 0; |
322 | } |
323 | if (!nSlacks) slackStart = -1; |
324 | return slackStart; |
325 | } |
326 | void |
327 | Idiot::crash(int numberPass, CoinMessageHandler * handler, |
328 | const CoinMessages *messages, bool doCrossover) |
329 | { |
330 | // lightweight options |
331 | int numberColumns = model_->getNumCols(); |
332 | const double * objective = model_->getObjCoefficients(); |
333 | int nnzero = 0; |
334 | double sum = 0.0; |
335 | int i; |
336 | for (i = 0; i < numberColumns; i++) { |
337 | if (objective[i]) { |
338 | sum += fabs(objective[i]); |
339 | nnzero++; |
340 | } |
341 | } |
342 | sum /= static_cast<double> (nnzero + 1); |
343 | if (maxIts_ == 5) |
344 | maxIts_ = 2; |
345 | if (numberPass <= 0) |
346 | majorIterations_ = static_cast<int>(2 + log10(static_cast<double>(numberColumns + 1))); |
347 | else |
348 | majorIterations_ = numberPass; |
349 | // If mu not changed then compute |
350 | if (mu_ == 1e-4) |
351 | mu_ = CoinMax(1.0e-3, sum * 1.0e-5); |
352 | if (maxIts2_ == 100) { |
353 | if (!lightWeight_) { |
354 | maxIts2_ = 105; |
355 | } else if (lightWeight_ == 1) { |
356 | mu_ *= 1000.0; |
357 | maxIts2_ = 23; |
358 | } else if (lightWeight_ == 2) { |
359 | maxIts2_ = 11; |
360 | } else { |
361 | maxIts2_ = 23; |
362 | } |
363 | } |
364 | //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_); |
365 | solve2(handler, messages); |
366 | #ifndef OSI_IDIOT |
367 | if (doCrossover) { |
368 | double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows()); |
369 | if ((averageInfeas < 0.01 && (strategy_ & 512) != 0) || (strategy_ & 8192) != 0) |
370 | crossOver(16 + 1); |
371 | else |
372 | crossOver(majorIterations_ < 1000000 ? 3 : 2); |
373 | } |
374 | #endif |
375 | } |
376 | |
377 | void |
378 | Idiot::solve() |
379 | { |
380 | CoinMessages dummy; |
381 | solve2(NULL, &dummy); |
382 | } |
383 | void |
384 | Idiot::solve2(CoinMessageHandler * handler, const CoinMessages * messages) |
385 | { |
386 | int strategy = 0; |
387 | double d2; |
388 | int i, n; |
389 | int allOnes = 1; |
390 | int iteration = 0; |
391 | int iterationTotal = 0; |
392 | int nTry = 0; /* number of tries at same weight */ |
393 | double fixTolerance = IDIOT_FIX_TOLERANCE; |
394 | int maxBigIts = maxBigIts_; |
395 | int maxIts = maxIts_; |
396 | int logLevel = logLevel_; |
397 | int saveMajorIterations = majorIterations_; |
398 | majorIterations_ = majorIterations_ % 1000000; |
399 | if (handler) { |
400 | if (handler->logLevel() > 0 && handler->logLevel() < 3) |
401 | logLevel = 1; |
402 | else if (!handler->logLevel()) |
403 | logLevel = 0; |
404 | else |
405 | logLevel = 7; |
406 | } |
407 | double djExit = djTolerance_; |
408 | double djFlag = 1.0 + 100.0 * djExit; |
409 | double djTol = 0.00001; |
410 | double mu = mu_; |
411 | double drop = drop_; |
412 | int maxIts2 = maxIts2_; |
413 | double factor = muFactor_; |
414 | double smallInfeas = smallInfeas_; |
415 | double reasonableInfeas = reasonableInfeas_; |
416 | double stopMu = stopMu_; |
417 | double maxmin, offset; |
418 | double lastWeighted = 1.0e50; |
419 | double exitDrop = exitDrop_; |
420 | double fakeSmall = smallInfeas; |
421 | double firstInfeas; |
422 | int badIts = 0; |
423 | int slackStart, slackEnd, ordStart, ordEnd; |
424 | int checkIteration = 0; |
425 | int lambdaIteration = 0; |
426 | int belowReasonable = 0; /* set if ever gone below reasonable infeas */ |
427 | double bestWeighted = 1.0e60; |
428 | double bestFeasible = 1.0e60; /* best solution while feasible */ |
429 | IdiotResult result, lastResult; |
430 | int saveStrategy = strategy_; |
431 | const int strategies[] = {0, 2, 128}; |
432 | double saveLambdaScale = 0.0; |
433 | if ((saveStrategy & 128) != 0) { |
434 | fixTolerance = SMALL_IDIOT_FIX_TOLERANCE; |
435 | } |
436 | #ifdef OSI_IDIOT |
437 | const CoinPackedMatrix * matrix = model_->getMatrixByCol(); |
438 | #else |
439 | ClpMatrixBase * matrix = model_->clpMatrix(); |
440 | #endif |
441 | const int * row = matrix->getIndices(); |
442 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
443 | const int * columnLength = matrix->getVectorLengths(); |
444 | const double * element = matrix->getElements(); |
445 | int nrows = model_->getNumRows(); |
446 | int ncols = model_->getNumCols(); |
447 | double * rowsol, * colsol; |
448 | double * pi, * dj; |
449 | #ifndef OSI_IDIOT |
450 | double * cost = model_->objective(); |
451 | double * lower = model_->columnLower(); |
452 | double * upper = model_->columnUpper(); |
453 | #else |
454 | double * cost = new double [ncols]; |
455 | CoinMemcpyN( model_->getObjCoefficients(), ncols, cost); |
456 | const double * lower = model_->getColLower(); |
457 | const double * upper = model_->getColUpper(); |
458 | #endif |
459 | const double *elemXX; |
460 | double * saveSol; |
461 | double * rowupper = new double[nrows]; // not const as modified |
462 | CoinMemcpyN(model_->getRowUpper(), nrows, rowupper); |
463 | double * rowlower = new double[nrows]; // not const as modified |
464 | CoinMemcpyN(model_->getRowLower(), nrows, rowlower); |
465 | CoinThreadRandom * randomNumberGenerator = model_->randomNumberGenerator(); |
466 | int * whenUsed; |
467 | double * lambda; |
468 | saveSol = new double[ncols]; |
469 | lambda = new double [nrows]; |
470 | rowsol = new double[nrows]; |
471 | colsol = new double [ncols]; |
472 | CoinMemcpyN(model_->getColSolution(), ncols, colsol); |
473 | pi = new double[nrows]; |
474 | dj = new double[ncols]; |
475 | delete [] whenUsed_; |
476 | bool oddSlacks = false; |
477 | // See if any costed slacks |
478 | int numberSlacks = 0; |
479 | for (i = 0; i < ncols; i++) { |
480 | if (columnLength[i] == 1) |
481 | numberSlacks++; |
482 | } |
483 | if (!numberSlacks) { |
484 | whenUsed_ = new int[ncols]; |
485 | } else { |
486 | #ifdef COIN_DEVELOP |
487 | printf("%d slacks\n" , numberSlacks); |
488 | #endif |
489 | oddSlacks = true; |
490 | int = static_cast<int> (nrows * sizeof(double) / sizeof(int)); |
491 | whenUsed_ = new int[2*ncols+2*nrows+extra]; |
492 | int * posSlack = whenUsed_ + ncols; |
493 | int * negSlack = posSlack + nrows; |
494 | int * nextSlack = negSlack + nrows; |
495 | for (i = 0; i < nrows; i++) { |
496 | posSlack[i] = -1; |
497 | negSlack[i] = -1; |
498 | } |
499 | for (i = 0; i < ncols; i++) |
500 | nextSlack[i] = -1; |
501 | for (i = 0; i < ncols; i++) { |
502 | if (columnLength[i] == 1) { |
503 | CoinBigIndex j = columnStart[i]; |
504 | int iRow = row[j]; |
505 | if (element[j] > 0.0) { |
506 | if (posSlack[iRow] == -1) { |
507 | posSlack[iRow] = i; |
508 | } else { |
509 | int iCol = posSlack[iRow]; |
510 | while (nextSlack[iCol] >= 0) |
511 | iCol = nextSlack[iCol]; |
512 | nextSlack[iCol] = i; |
513 | } |
514 | } else { |
515 | if (negSlack[iRow] == -1) { |
516 | negSlack[iRow] = i; |
517 | } else { |
518 | int iCol = negSlack[iRow]; |
519 | while (nextSlack[iCol] >= 0) |
520 | iCol = nextSlack[iCol]; |
521 | nextSlack[iCol] = i; |
522 | } |
523 | } |
524 | } |
525 | } |
526 | // now sort |
527 | for (i = 0; i < nrows; i++) { |
528 | int iCol; |
529 | iCol = posSlack[i]; |
530 | if (iCol >= 0) { |
531 | CoinBigIndex j = columnStart[iCol]; |
532 | #ifndef NDEBUG |
533 | int iRow = row[j]; |
534 | #endif |
535 | assert (element[j] > 0.0); |
536 | assert (iRow == i); |
537 | dj[0] = cost[iCol] / element[j]; |
538 | whenUsed_[0] = iCol; |
539 | int n = 1; |
540 | while (nextSlack[iCol] >= 0) { |
541 | iCol = nextSlack[iCol]; |
542 | CoinBigIndex j = columnStart[iCol]; |
543 | #ifndef NDEBUG |
544 | int iRow = row[j]; |
545 | #endif |
546 | assert (element[j] > 0.0); |
547 | assert (iRow == i); |
548 | dj[n] = cost[iCol] / element[j]; |
549 | whenUsed_[n++] = iCol; |
550 | } |
551 | for (j = 0; j < n; j++) { |
552 | int jCol = whenUsed_[j]; |
553 | nextSlack[jCol] = -2; |
554 | } |
555 | CoinSort_2(dj, dj + n, whenUsed_); |
556 | // put back |
557 | iCol = whenUsed_[0]; |
558 | posSlack[i] = iCol; |
559 | for (j = 1; j < n; j++) { |
560 | int jCol = whenUsed_[j]; |
561 | nextSlack[iCol] = jCol; |
562 | iCol = jCol; |
563 | } |
564 | } |
565 | iCol = negSlack[i]; |
566 | if (iCol >= 0) { |
567 | CoinBigIndex j = columnStart[iCol]; |
568 | #ifndef NDEBUG |
569 | int iRow = row[j]; |
570 | #endif |
571 | assert (element[j] < 0.0); |
572 | assert (iRow == i); |
573 | dj[0] = -cost[iCol] / element[j]; |
574 | whenUsed_[0] = iCol; |
575 | int n = 1; |
576 | while (nextSlack[iCol] >= 0) { |
577 | iCol = nextSlack[iCol]; |
578 | CoinBigIndex j = columnStart[iCol]; |
579 | #ifndef NDEBUG |
580 | int iRow = row[j]; |
581 | #endif |
582 | assert (element[j] < 0.0); |
583 | assert (iRow == i); |
584 | dj[n] = -cost[iCol] / element[j]; |
585 | whenUsed_[n++] = iCol; |
586 | } |
587 | for (j = 0; j < n; j++) { |
588 | int jCol = whenUsed_[j]; |
589 | nextSlack[jCol] = -2; |
590 | } |
591 | CoinSort_2(dj, dj + n, whenUsed_); |
592 | // put back |
593 | iCol = whenUsed_[0]; |
594 | negSlack[i] = iCol; |
595 | for (j = 1; j < n; j++) { |
596 | int jCol = whenUsed_[j]; |
597 | nextSlack[iCol] = jCol; |
598 | iCol = jCol; |
599 | } |
600 | } |
601 | } |
602 | } |
603 | whenUsed = whenUsed_; |
604 | if (model_->getObjSense() == -1.0) { |
605 | maxmin = -1.0; |
606 | } else { |
607 | maxmin = 1.0; |
608 | } |
609 | model_->getDblParam(OsiObjOffset, offset); |
610 | if (!maxIts2) maxIts2 = maxIts; |
611 | strategy = strategy_; |
612 | strategy &= 3; |
613 | memset(lambda, 0, nrows * sizeof(double)); |
614 | slackStart = countCostedSlacks(model_); |
615 | if (slackStart >= 0) { |
616 | COIN_DETAIL_PRINT(printf("This model has costed slacks\n" )); |
617 | slackEnd = slackStart + nrows; |
618 | if (slackStart) { |
619 | ordStart = 0; |
620 | ordEnd = slackStart; |
621 | } else { |
622 | ordStart = nrows; |
623 | ordEnd = ncols; |
624 | } |
625 | } else { |
626 | slackEnd = slackStart; |
627 | ordStart = 0; |
628 | ordEnd = ncols; |
629 | } |
630 | if (offset && logLevel > 2) { |
631 | printf("** Objective offset is %g\n" , offset); |
632 | } |
633 | /* compute reasonable solution cost */ |
634 | for (i = 0; i < nrows; i++) { |
635 | rowsol[i] = 1.0e31; |
636 | } |
637 | for (i = 0; i < ncols; i++) { |
638 | CoinBigIndex j; |
639 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
640 | if (element[j] != 1.0) { |
641 | allOnes = 0; |
642 | break; |
643 | } |
644 | } |
645 | } |
646 | if (allOnes) { |
647 | elemXX = NULL; |
648 | } else { |
649 | elemXX = element; |
650 | } |
651 | // Do scaling if wanted |
652 | bool scaled = false; |
653 | #ifndef OSI_IDIOT |
654 | if ((strategy_ & 32) != 0 && !allOnes) { |
655 | if (model_->scalingFlag() > 0) |
656 | scaled = model_->clpMatrix()->scale(model_) == 0; |
657 | if (scaled) { |
658 | const double * rowScale = model_->rowScale(); |
659 | const double * columnScale = model_->columnScale(); |
660 | double * oldLower = lower; |
661 | double * oldUpper = upper; |
662 | double * oldCost = cost; |
663 | lower = new double[ncols]; |
664 | upper = new double[ncols]; |
665 | cost = new double[ncols]; |
666 | CoinMemcpyN(oldLower, ncols, lower); |
667 | CoinMemcpyN(oldUpper, ncols, upper); |
668 | CoinMemcpyN(oldCost, ncols, cost); |
669 | int icol, irow; |
670 | for (icol = 0; icol < ncols; icol++) { |
671 | double multiplier = 1.0 / columnScale[icol]; |
672 | if (lower[icol] > -1.0e50) |
673 | lower[icol] *= multiplier; |
674 | if (upper[icol] < 1.0e50) |
675 | upper[icol] *= multiplier; |
676 | colsol[icol] *= multiplier; |
677 | cost[icol] *= columnScale[icol]; |
678 | } |
679 | CoinMemcpyN(model_->rowLower(), nrows, rowlower); |
680 | for (irow = 0; irow < nrows; irow++) { |
681 | double multiplier = rowScale[irow]; |
682 | if (rowlower[irow] > -1.0e50) |
683 | rowlower[irow] *= multiplier; |
684 | if (rowupper[irow] < 1.0e50) |
685 | rowupper[irow] *= multiplier; |
686 | rowsol[irow] *= multiplier; |
687 | } |
688 | int length = columnStart[ncols-1] + columnLength[ncols-1]; |
689 | double * elemYY = new double[length]; |
690 | for (i = 0; i < ncols; i++) { |
691 | CoinBigIndex j; |
692 | double scale = columnScale[i]; |
693 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
694 | int irow = row[j]; |
695 | elemYY[j] = element[j] * scale * rowScale[irow]; |
696 | } |
697 | } |
698 | elemXX = elemYY; |
699 | } |
700 | } |
701 | #endif |
702 | for (i = 0; i < ncols; i++) { |
703 | CoinBigIndex j; |
704 | double dd = columnLength[i]; |
705 | dd = cost[i] / dd; |
706 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
707 | int irow = row[j]; |
708 | if (dd < rowsol[irow]) { |
709 | rowsol[irow] = dd; |
710 | } |
711 | } |
712 | } |
713 | d2 = 0.0; |
714 | for (i = 0; i < nrows; i++) { |
715 | d2 += rowsol[i]; |
716 | } |
717 | d2 *= 2.0; /* for luck */ |
718 | |
719 | d2 = d2 / static_cast<double> (4 * nrows + 8000); |
720 | d2 *= 0.5; /* halve with more flexible method */ |
721 | if (d2 < 5.0) d2 = 5.0; |
722 | if (djExit == 0.0) { |
723 | djExit = d2; |
724 | } |
725 | if ((saveStrategy & 4) != 0) { |
726 | /* go to relative tolerances - first small */ |
727 | djExit = 1.0e-10; |
728 | djFlag = 1.0e-5; |
729 | drop = 1.0e-10; |
730 | } |
731 | memset(whenUsed, 0, ncols * sizeof(int)); |
732 | strategy = strategies[strategy]; |
733 | if ((saveStrategy & 8) != 0) strategy |= 64; /* don't allow large theta */ |
734 | CoinMemcpyN(colsol, ncols, saveSol); |
735 | |
736 | lastResult = IdiSolve(nrows, ncols, rowsol , colsol, pi, |
737 | dj, cost, rowlower, rowupper, |
738 | lower, upper, elemXX, row, columnStart, columnLength, lambda, |
739 | 0, mu, drop, |
740 | maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator); |
741 | // update whenUsed_ |
742 | n = cleanIteration(iteration, ordStart, ordEnd, |
743 | colsol, lower, upper, |
744 | rowlower, rowupper, |
745 | cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas); |
746 | if ((strategy_ & 16384) != 0) { |
747 | int * posSlack = whenUsed_ + ncols; |
748 | int * negSlack = posSlack + nrows; |
749 | int * nextSlack = negSlack + nrows; |
750 | double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols); |
751 | for (i = 0; i < nrows; i++) |
752 | rowsol[i] += rowsol2[i]; |
753 | } |
754 | if ((logLevel_ & 1) != 0) { |
755 | #ifndef OSI_IDIOT |
756 | if (!handler) { |
757 | #endif |
758 | printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n" , |
759 | iteration, lastResult.infeas, lastResult.objval, mu, lastResult.iteration, n); |
760 | #ifndef OSI_IDIOT |
761 | } else { |
762 | handler->message(CLP_IDIOT_ITERATION, *messages) |
763 | << iteration << lastResult.infeas << lastResult.objval << mu << lastResult.iteration << n |
764 | << CoinMessageEol; |
765 | } |
766 | #endif |
767 | } |
768 | int numberBaseTrys = 0; // for first time |
769 | int numberAway = -1; |
770 | iterationTotal = lastResult.iteration; |
771 | firstInfeas = lastResult.infeas; |
772 | if ((strategy_ & 1024) != 0) reasonableInfeas = 0.5 * firstInfeas; |
773 | if (lastResult.infeas < reasonableInfeas) lastResult.infeas = reasonableInfeas; |
774 | double keepinfeas = 1.0e31; |
775 | double lastInfeas = 1.0e31; |
776 | double bestInfeas = 1.0e31; |
777 | while ((mu > stopMu && lastResult.infeas > smallInfeas) || |
778 | (lastResult.infeas <= smallInfeas && |
779 | dropping(lastResult, exitDrop, smallInfeas, &badIts)) || |
780 | checkIteration < 2 || lambdaIteration < lambdaIterations_) { |
781 | if (lastResult.infeas <= exitFeasibility_) |
782 | break; |
783 | iteration++; |
784 | checkIteration++; |
785 | if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) { |
786 | bestFeasible = lastResult.objval; |
787 | } |
788 | if (lastResult.infeas + mu * lastResult.objval < bestWeighted) { |
789 | bestWeighted = lastResult.objval + mu * lastResult.objval; |
790 | } |
791 | if ((saveStrategy & 4096)) strategy |= 256; |
792 | if ((saveStrategy & 4) != 0 && iteration > 2) { |
793 | /* go to relative tolerances */ |
794 | double weighted = 10.0 + fabs(lastWeighted); |
795 | djExit = djTolerance_ * weighted; |
796 | djFlag = 2.0 * djExit; |
797 | drop = drop_ * weighted; |
798 | djTol = 0.01 * djExit; |
799 | } |
800 | CoinMemcpyN(colsol, ncols, saveSol); |
801 | result = IdiSolve(nrows, ncols, rowsol , colsol, pi, dj, |
802 | cost, rowlower, rowupper, |
803 | lower, upper, elemXX, row, columnStart, columnLength, lambda, |
804 | maxIts, mu, drop, |
805 | maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator); |
806 | n = cleanIteration(iteration, ordStart, ordEnd, |
807 | colsol, lower, upper, |
808 | rowlower, rowupper, |
809 | cost, elemXX, fixTolerance, result.objval, result.infeas); |
810 | if ((strategy_ & 16384) != 0) { |
811 | int * posSlack = whenUsed_ + ncols; |
812 | int * negSlack = posSlack + nrows; |
813 | int * nextSlack = negSlack + nrows; |
814 | double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols); |
815 | for (i = 0; i < nrows; i++) |
816 | rowsol[i] += rowsol2[i]; |
817 | } |
818 | if ((logLevel_ & 1) != 0) { |
819 | #ifndef OSI_IDIOT |
820 | if (!handler) { |
821 | #endif |
822 | printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n" , |
823 | iteration, result.infeas, result.objval, mu, result.iteration, n); |
824 | #ifndef OSI_IDIOT |
825 | } else { |
826 | handler->message(CLP_IDIOT_ITERATION, *messages) |
827 | << iteration << result.infeas << result.objval << mu << result.iteration << n |
828 | << CoinMessageEol; |
829 | } |
830 | #endif |
831 | } |
832 | if (iteration > 50 && n == numberAway && result.infeas < 1.0e-4) { |
833 | #ifdef CLP_INVESTIGATE |
834 | printf("infeas small %g\n" , result.infeas); |
835 | #endif |
836 | break; // not much happening |
837 | } |
838 | if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) { |
839 | if (lastInfeas != bestInfeas && CoinMin(result.infeas, lastInfeas) > 0.95 * bestInfeas) |
840 | majorIterations_ = CoinMin(majorIterations_, iteration); // not getting feasible |
841 | } |
842 | lastInfeas = result.infeas; |
843 | numberAway = n; |
844 | keepinfeas = result.infeas; |
845 | lastWeighted = result.weighted; |
846 | iterationTotal += result.iteration; |
847 | if (iteration == 1) { |
848 | if ((strategy_ & 1024) != 0 && mu < 1.0e-10) |
849 | result.infeas = firstInfeas * 0.8; |
850 | if (majorIterations_ >= 50 || dropEnoughFeasibility_ <= 0.0) |
851 | result.infeas *= 0.8; |
852 | if (result.infeas > firstInfeas * 0.9 |
853 | && result.infeas > reasonableInfeas) { |
854 | iteration--; |
855 | if (majorIterations_ < 50) |
856 | mu *= 1.0e-1; |
857 | else |
858 | mu *= 0.7; |
859 | bestFeasible = 1.0e31; |
860 | bestWeighted = 1.0e60; |
861 | numberBaseTrys++; |
862 | if (mu < 1.0e-30 || (numberBaseTrys > 10 && lightWeight_)) { |
863 | // back to all slack basis |
864 | lightWeight_ = 2; |
865 | break; |
866 | } |
867 | CoinMemcpyN(saveSol, ncols, colsol); |
868 | } else { |
869 | maxIts = maxIts2; |
870 | checkIteration = 0; |
871 | if ((strategy_ & 1024) != 0) mu *= 1.0e-1; |
872 | } |
873 | } else { |
874 | } |
875 | bestInfeas = CoinMin(bestInfeas, result.infeas); |
876 | if (iteration) { |
877 | /* this code is in to force it to terminate sometime */ |
878 | double changeMu = factor; |
879 | if ((saveStrategy & 64) != 0) { |
880 | keepinfeas = 0.0; /* switch off ranga's increase */ |
881 | fakeSmall = smallInfeas; |
882 | } else { |
883 | fakeSmall = -1.0; |
884 | } |
885 | saveLambdaScale = 0.0; |
886 | if (result.infeas > reasonableInfeas || |
887 | (nTry + 1 == maxBigIts && result.infeas > fakeSmall)) { |
888 | if (result.infeas > lastResult.infeas*(1.0 - dropEnoughFeasibility_) || |
889 | nTry + 1 == maxBigIts || |
890 | (result.infeas > lastResult.infeas * 0.9 |
891 | && result.weighted > lastResult.weighted |
892 | - dropEnoughWeighted_ * CoinMax(fabs(lastResult.weighted), fabs(result.weighted)))) { |
893 | mu *= changeMu; |
894 | if ((saveStrategy & 32) != 0 && result.infeas < reasonableInfeas && 0) { |
895 | reasonableInfeas = CoinMax(smallInfeas, reasonableInfeas * sqrt(changeMu)); |
896 | COIN_DETAIL_PRINT(printf("reasonable infeas now %g\n" , reasonableInfeas)); |
897 | } |
898 | result.weighted = 1.0e60; |
899 | nTry = 0; |
900 | bestFeasible = 1.0e31; |
901 | bestWeighted = 1.0e60; |
902 | checkIteration = 0; |
903 | lambdaIteration = 0; |
904 | #define LAMBDA |
905 | #ifdef LAMBDA |
906 | if ((saveStrategy & 2048) == 0) { |
907 | memset(lambda, 0, nrows * sizeof(double)); |
908 | } |
909 | #else |
910 | memset(lambda, 0, nrows * sizeof(double)); |
911 | #endif |
912 | } else { |
913 | nTry++; |
914 | } |
915 | } else if (lambdaIterations_ >= 0) { |
916 | /* update lambda */ |
917 | double scale = 1.0 / mu; |
918 | int i, nnz = 0; |
919 | saveLambdaScale = scale; |
920 | lambdaIteration++; |
921 | if ((saveStrategy & 4) == 0) drop = drop_ / 50.0; |
922 | if (lambdaIteration > 4 && |
923 | (((lambdaIteration % 10) == 0 && smallInfeas < keepinfeas) || |
924 | ((lambdaIteration % 5) == 0 && 1.5 * smallInfeas < keepinfeas))) { |
925 | //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas); |
926 | smallInfeas *= 1.5; |
927 | } |
928 | if ((saveStrategy & 2048) == 0) { |
929 | for (i = 0; i < nrows; i++) { |
930 | if (lambda[i]) nnz++; |
931 | lambda[i] += scale * rowsol[i]; |
932 | } |
933 | } else { |
934 | nnz = 1; |
935 | #ifdef LAMBDA |
936 | for (i = 0; i < nrows; i++) { |
937 | lambda[i] += scale * rowsol[i]; |
938 | } |
939 | #else |
940 | for (i = 0; i < nrows; i++) { |
941 | lambda[i] = scale * rowsol[i]; |
942 | } |
943 | for (i = 0; i < ncols; i++) { |
944 | CoinBigIndex j; |
945 | double value = cost[i] * maxmin; |
946 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
947 | int irow = row[j]; |
948 | value += element[j] * lambda[irow]; |
949 | } |
950 | cost[i] = value * maxmin; |
951 | } |
952 | for (i = 0; i < nrows; i++) { |
953 | offset += lambda[i] * rowupper[i]; |
954 | lambda[i] = 0.0; |
955 | } |
956 | #ifdef DEBUG |
957 | printf("offset %g\n" , offset); |
958 | #endif |
959 | model_->setDblParam(OsiObjOffset, offset); |
960 | #endif |
961 | } |
962 | nTry++; |
963 | if (!nnz) { |
964 | bestFeasible = 1.0e32; |
965 | bestWeighted = 1.0e60; |
966 | checkIteration = 0; |
967 | result.weighted = 1.0e31; |
968 | } |
969 | #ifdef DEBUG |
970 | double trueCost = 0.0; |
971 | for (i = 0; i < ncols; i++) { |
972 | int j; |
973 | trueCost += cost[i] * colsol[i]; |
974 | } |
975 | printf("True objective %g\n" , trueCost - offset); |
976 | #endif |
977 | } else { |
978 | nTry++; |
979 | } |
980 | lastResult = result; |
981 | if (result.infeas < reasonableInfeas && !belowReasonable) { |
982 | belowReasonable = 1; |
983 | bestFeasible = 1.0e32; |
984 | bestWeighted = 1.0e60; |
985 | checkIteration = 0; |
986 | result.weighted = 1.0e31; |
987 | } |
988 | } |
989 | if (iteration >= majorIterations_) { |
990 | // If not feasible and crash then dive dive dive |
991 | if (mu > 1.0e-12 && result.infeas > 1.0 && majorIterations_ < 40) { |
992 | mu = 1.0e-30; |
993 | majorIterations_ = iteration + 1; |
994 | stopMu = 0.0; |
995 | } else { |
996 | if (logLevel > 2) |
997 | printf("Exiting due to number of major iterations\n" ); |
998 | break; |
999 | } |
1000 | } |
1001 | } |
1002 | majorIterations_ = saveMajorIterations; |
1003 | #ifndef OSI_IDIOT |
1004 | if (scaled) { |
1005 | // Scale solution and free arrays |
1006 | const double * rowScale = model_->rowScale(); |
1007 | const double * columnScale = model_->columnScale(); |
1008 | int icol, irow; |
1009 | for (icol = 0; icol < ncols; icol++) { |
1010 | colsol[icol] *= columnScale[icol]; |
1011 | saveSol[icol] *= columnScale[icol]; |
1012 | dj[icol] /= columnScale[icol]; |
1013 | } |
1014 | for (irow = 0; irow < nrows; irow++) { |
1015 | rowsol[irow] /= rowScale[irow]; |
1016 | pi[irow] *= rowScale[irow]; |
1017 | } |
1018 | // Don't know why getting Microsoft problems |
1019 | #if defined (_MSC_VER) |
1020 | delete [] ( double *) elemXX; |
1021 | #else |
1022 | delete [] elemXX; |
1023 | #endif |
1024 | model_->setRowScale(NULL); |
1025 | model_->setColumnScale(NULL); |
1026 | delete [] lower; |
1027 | delete [] upper; |
1028 | delete [] cost; |
1029 | lower = model_->columnLower(); |
1030 | upper = model_->columnUpper(); |
1031 | cost = model_->objective(); |
1032 | //rowlower = model_->rowLower(); |
1033 | } |
1034 | #endif |
1035 | #define TRYTHIS |
1036 | #ifdef TRYTHIS |
1037 | if ((saveStrategy & 2048) != 0) { |
1038 | double offset; |
1039 | model_->getDblParam(OsiObjOffset, offset); |
1040 | for (i = 0; i < ncols; i++) { |
1041 | CoinBigIndex j; |
1042 | double djval = cost[i] * maxmin; |
1043 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
1044 | int irow = row[j]; |
1045 | djval -= element[j] * lambda[irow]; |
1046 | } |
1047 | cost[i] = djval; |
1048 | } |
1049 | for (i = 0; i < nrows; i++) { |
1050 | offset += lambda[i] * rowupper[i]; |
1051 | } |
1052 | model_->setDblParam(OsiObjOffset, offset); |
1053 | } |
1054 | #endif |
1055 | if (saveLambdaScale) { |
1056 | /* back off last update */ |
1057 | for (i = 0; i < nrows; i++) { |
1058 | lambda[i] -= saveLambdaScale * rowsol[i]; |
1059 | } |
1060 | } |
1061 | muAtExit_ = mu; |
1062 | // For last iteration make as feasible as possible |
1063 | if (oddSlacks) |
1064 | strategy_ |= 16384; |
1065 | // not scaled |
1066 | n = cleanIteration(iteration, ordStart, ordEnd, |
1067 | colsol, lower, upper, |
1068 | model_->rowLower(), model_->rowUpper(), |
1069 | cost, element, fixTolerance, lastResult.objval, lastResult.infeas); |
1070 | #if 0 |
1071 | if ((logLevel & 1) == 0 || (strategy_ & 16384) != 0) { |
1072 | printf( |
1073 | "%d - mu %g, infeasibility %g, objective %g, %d interior\n" , |
1074 | iteration, mu, lastResult.infeas, lastResult.objval, n); |
1075 | } |
1076 | #endif |
1077 | #ifndef OSI_IDIOT |
1078 | model_->setSumPrimalInfeasibilities(lastResult.infeas); |
1079 | #endif |
1080 | // Put back more feasible solution |
1081 | double saveInfeas[] = {0.0, 0.0}; |
1082 | for (int iSol = 0; iSol < 3; iSol++) { |
1083 | const double * solution = iSol ? colsol : saveSol; |
1084 | if (iSol == 2 && saveInfeas[0] < saveInfeas[1]) { |
1085 | // put back best solution |
1086 | CoinMemcpyN(saveSol, ncols, colsol); |
1087 | } |
1088 | double large = 0.0; |
1089 | int i; |
1090 | memset(rowsol, 0, nrows * sizeof(double)); |
1091 | for (i = 0; i < ncols; i++) { |
1092 | CoinBigIndex j; |
1093 | double value = solution[i]; |
1094 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
1095 | int irow = row[j]; |
1096 | rowsol[irow] += element[j] * value; |
1097 | } |
1098 | } |
1099 | for (i = 0; i < nrows; i++) { |
1100 | if (rowsol[i] > rowupper[i]) { |
1101 | double diff = rowsol[i] - rowupper[i]; |
1102 | if (diff > large) |
1103 | large = diff; |
1104 | } else if (rowsol[i] < rowlower[i]) { |
1105 | double diff = rowlower[i] - rowsol[i]; |
1106 | if (diff > large) |
1107 | large = diff; |
1108 | } |
1109 | } |
1110 | if (iSol < 2) |
1111 | saveInfeas[iSol] = large; |
1112 | if (logLevel > 2) |
1113 | printf("largest infeasibility is %g\n" , large); |
1114 | } |
1115 | /* subtract out lambda */ |
1116 | for (i = 0; i < nrows; i++) { |
1117 | pi[i] -= lambda[i]; |
1118 | } |
1119 | for (i = 0; i < ncols; i++) { |
1120 | CoinBigIndex j; |
1121 | double djval = cost[i] * maxmin; |
1122 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
1123 | int irow = row[j]; |
1124 | djval -= element[j] * pi[irow]; |
1125 | } |
1126 | dj[i] = djval; |
1127 | } |
1128 | if ((strategy_ & 1024) != 0) { |
1129 | double ratio = static_cast<double> (ncols) / static_cast<double> (nrows); |
1130 | COIN_DETAIL_PRINT(printf("col/row ratio %g infeas ratio %g\n" , ratio, lastResult.infeas / firstInfeas)); |
1131 | if (lastResult.infeas > 0.01 * firstInfeas * ratio) { |
1132 | strategy_ &= (~1024); |
1133 | COIN_DETAIL_PRINT(printf(" - layer off\n" )); |
1134 | } else { |
1135 | COIN_DETAIL_PRINT(printf(" - layer on\n" )); |
1136 | } |
1137 | } |
1138 | delete [] saveSol; |
1139 | delete [] lambda; |
1140 | // save solution |
1141 | // duals not much use - but save anyway |
1142 | #ifndef OSI_IDIOT |
1143 | CoinMemcpyN(rowsol, nrows, model_->primalRowSolution()); |
1144 | CoinMemcpyN(colsol, ncols, model_->primalColumnSolution()); |
1145 | CoinMemcpyN(pi, nrows, model_->dualRowSolution()); |
1146 | CoinMemcpyN(dj, ncols, model_->dualColumnSolution()); |
1147 | #else |
1148 | model_->setColSolution(colsol); |
1149 | model_->setRowPrice(pi); |
1150 | delete [] cost; |
1151 | #endif |
1152 | delete [] rowsol; |
1153 | delete [] colsol; |
1154 | delete [] pi; |
1155 | delete [] dj; |
1156 | delete [] rowlower; |
1157 | delete [] rowupper; |
1158 | return ; |
1159 | } |
1160 | #ifndef OSI_IDIOT |
1161 | void |
1162 | Idiot::crossOver(int mode) |
1163 | { |
1164 | if (lightWeight_ == 2) { |
1165 | // total failure |
1166 | model_->allSlackBasis(); |
1167 | return; |
1168 | } |
1169 | double fixTolerance = IDIOT_FIX_TOLERANCE; |
1170 | #ifdef COIN_DEVELOP |
1171 | double startTime = CoinCpuTime(); |
1172 | #endif |
1173 | ClpSimplex * saveModel = NULL; |
1174 | ClpMatrixBase * matrix = model_->clpMatrix(); |
1175 | const int * row = matrix->getIndices(); |
1176 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
1177 | const int * columnLength = matrix->getVectorLengths(); |
1178 | const double * element = matrix->getElements(); |
1179 | const double * rowupper = model_->getRowUpper(); |
1180 | int nrows = model_->getNumRows(); |
1181 | int ncols = model_->getNumCols(); |
1182 | double * rowsol, * colsol; |
1183 | // different for Osi |
1184 | double * lower = model_->columnLower(); |
1185 | double * upper = model_->columnUpper(); |
1186 | const double * rowlower = model_->getRowLower(); |
1187 | int * whenUsed = whenUsed_; |
1188 | rowsol = model_->primalRowSolution(); |
1189 | colsol = model_->primalColumnSolution(); |
1190 | double * cost = model_->objective(); |
1191 | |
1192 | int slackEnd, ordStart, ordEnd; |
1193 | int slackStart = countCostedSlacks(model_); |
1194 | |
1195 | int addAll = mode & 7; |
1196 | int presolve = 0; |
1197 | |
1198 | double djTolerance = djTolerance_; |
1199 | if (djTolerance > 0.0 && djTolerance < 1.0) |
1200 | djTolerance = 1.0; |
1201 | int iteration; |
1202 | int i, n = 0; |
1203 | double ratio = 1.0; |
1204 | double objValue = 0.0; |
1205 | if ((strategy_ & 128) != 0) { |
1206 | fixTolerance = SMALL_IDIOT_FIX_TOLERANCE; |
1207 | } |
1208 | if ((mode & 16) != 0 && addAll < 3) presolve = 1; |
1209 | double * saveUpper = NULL; |
1210 | double * saveLower = NULL; |
1211 | double * saveRowUpper = NULL; |
1212 | double * saveRowLower = NULL; |
1213 | bool allowInfeasible = ((strategy_ & 8192) != 0) || (majorIterations_ > 1000000); |
1214 | if (addAll < 3) { |
1215 | saveUpper = new double [ncols]; |
1216 | saveLower = new double [ncols]; |
1217 | CoinMemcpyN(upper, ncols, saveUpper); |
1218 | CoinMemcpyN(lower, ncols, saveLower); |
1219 | if (allowInfeasible) { |
1220 | saveRowUpper = new double [nrows]; |
1221 | saveRowLower = new double [nrows]; |
1222 | CoinMemcpyN(rowupper, nrows, saveRowUpper); |
1223 | CoinMemcpyN(rowlower, nrows, saveRowLower); |
1224 | double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows()); |
1225 | fixTolerance = CoinMax(fixTolerance, 1.0e-5 * averageInfeas); |
1226 | } |
1227 | } |
1228 | if (slackStart >= 0) { |
1229 | slackEnd = slackStart + nrows; |
1230 | if (slackStart) { |
1231 | ordStart = 0; |
1232 | ordEnd = slackStart; |
1233 | } else { |
1234 | ordStart = nrows; |
1235 | ordEnd = ncols; |
1236 | } |
1237 | } else { |
1238 | slackEnd = slackStart; |
1239 | ordStart = 0; |
1240 | ordEnd = ncols; |
1241 | } |
1242 | /* get correct rowsol (without known slacks) */ |
1243 | memset(rowsol, 0, nrows * sizeof(double)); |
1244 | for (i = ordStart; i < ordEnd; i++) { |
1245 | CoinBigIndex j; |
1246 | double value = colsol[i]; |
1247 | if (value < lower[i] + fixTolerance) { |
1248 | value = lower[i]; |
1249 | colsol[i] = value; |
1250 | } |
1251 | for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) { |
1252 | int irow = row[j]; |
1253 | rowsol[irow] += value * element[j]; |
1254 | } |
1255 | } |
1256 | if (slackStart >= 0) { |
1257 | for (i = 0; i < nrows; i++) { |
1258 | if (ratio * rowsol[i] > rowlower[i] && rowsol[i] > 1.0e-8) { |
1259 | ratio = rowlower[i] / rowsol[i]; |
1260 | } |
1261 | } |
1262 | for (i = 0; i < nrows; i++) { |
1263 | rowsol[i] *= ratio; |
1264 | } |
1265 | for (i = ordStart; i < ordEnd; i++) { |
1266 | double value = colsol[i] * ratio; |
1267 | colsol[i] = value; |
1268 | objValue += value * cost[i]; |
1269 | } |
1270 | for (i = 0; i < nrows; i++) { |
1271 | double value = rowlower[i] - rowsol[i]; |
1272 | colsol[i+slackStart] = value; |
1273 | objValue += value * cost[i+slackStart]; |
1274 | } |
1275 | COIN_DETAIL_PRINT(printf("New objective after scaling %g\n" , objValue)); |
1276 | } |
1277 | #if 0 |
1278 | maybe put back - but just get feasible ? |
1279 | // If not many fixed then just exit |
1280 | int numberFixed = 0; |
1281 | for (i = ordStart; i < ordEnd; i++) { |
1282 | if (colsol[i] < lower[i] + fixTolerance) |
1283 | numberFixed++; |
1284 | else if (colsol[i] > upper[i] - fixTolerance) |
1285 | numberFixed++; |
1286 | } |
1287 | if (numberFixed < ncols / 2) { |
1288 | addAll = 3; |
1289 | presolve = 0; |
1290 | } |
1291 | #endif |
1292 | #ifdef FEB_TRY |
1293 | int savePerturbation = model_->perturbation(); |
1294 | int saveOptions = model_->specialOptions(); |
1295 | model_->setSpecialOptions(saveOptions | 8192); |
1296 | if (savePerturbation_ == 50) |
1297 | model_->setPerturbation(56); |
1298 | #endif |
1299 | model_->createStatus(); |
1300 | /* addAll |
1301 | 0 - chosen,all used, all |
1302 | 1 - chosen, all |
1303 | 2 - all |
1304 | 3 - do not do anything - maybe basis |
1305 | */ |
1306 | for (i = ordStart; i < ordEnd; i++) { |
1307 | if (addAll < 2) { |
1308 | if (colsol[i] < lower[i] + fixTolerance) { |
1309 | upper[i] = lower[i]; |
1310 | colsol[i] = lower[i]; |
1311 | } else if (colsol[i] > upper[i] - fixTolerance) { |
1312 | lower[i] = upper[i]; |
1313 | colsol[i] = upper[i]; |
1314 | } |
1315 | } |
1316 | model_->setColumnStatus(i, ClpSimplex::superBasic); |
1317 | } |
1318 | if ((strategy_ & 16384) != 0) { |
1319 | // put in basis |
1320 | int * posSlack = whenUsed_ + ncols; |
1321 | int * negSlack = posSlack + nrows; |
1322 | int * nextSlack = negSlack + nrows; |
1323 | /* Laci - try both ways - to see what works - |
1324 | you can change second part as much as you want */ |
1325 | #ifndef LACI_TRY // was #if 1 |
1326 | // Array for sorting out slack values |
1327 | double * ratio = new double [ncols]; |
1328 | int * which = new int [ncols]; |
1329 | for (i = 0; i < nrows; i++) { |
1330 | if (posSlack[i] >= 0 || negSlack[i] >= 0) { |
1331 | int iCol; |
1332 | int nPlus = 0; |
1333 | int nMinus = 0; |
1334 | bool possible = true; |
1335 | // Get sum |
1336 | double sum = 0.0; |
1337 | iCol = posSlack[i]; |
1338 | while (iCol >= 0) { |
1339 | double value = element[columnStart[iCol]]; |
1340 | sum += value * colsol[iCol]; |
1341 | if (lower[iCol]) { |
1342 | possible = false; |
1343 | break; |
1344 | } else { |
1345 | nPlus++; |
1346 | } |
1347 | iCol = nextSlack[iCol]; |
1348 | } |
1349 | iCol = negSlack[i]; |
1350 | while (iCol >= 0) { |
1351 | double value = -element[columnStart[iCol]]; |
1352 | sum -= value * colsol[iCol]; |
1353 | if (lower[iCol]) { |
1354 | possible = false; |
1355 | break; |
1356 | } else { |
1357 | nMinus++; |
1358 | } |
1359 | iCol = nextSlack[iCol]; |
1360 | } |
1361 | //printf("%d plus, %d minus",nPlus,nMinus); |
1362 | //printf("\n"); |
1363 | if ((rowsol[i] - rowlower[i] < 1.0e-7 || |
1364 | rowupper[i] - rowsol[i] < 1.0e-7) && |
1365 | nPlus + nMinus < 2) |
1366 | possible = false; |
1367 | if (possible) { |
1368 | // Amount contributed by other varaibles |
1369 | sum = rowsol[i] - sum; |
1370 | double lo = rowlower[i]; |
1371 | if (lo > -1.0e20) |
1372 | lo -= sum; |
1373 | double up = rowupper[i]; |
1374 | if (up < 1.0e20) |
1375 | up -= sum; |
1376 | //printf("row bounds %g %g\n",lo,up); |
1377 | if (0) { |
1378 | double sum = 0.0; |
1379 | double x = 0.0; |
1380 | for (int k = 0; k < ncols; k++) { |
1381 | CoinBigIndex j; |
1382 | double value = colsol[k]; |
1383 | x += value * cost[k]; |
1384 | for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) { |
1385 | int irow = row[j]; |
1386 | if (irow == i) |
1387 | sum += element[j] * value; |
1388 | } |
1389 | } |
1390 | printf("Before sum %g <= %g <= %g cost %.18g\n" , |
1391 | rowlower[i], sum, rowupper[i], x); |
1392 | } |
1393 | // set all to zero |
1394 | iCol = posSlack[i]; |
1395 | while (iCol >= 0) { |
1396 | colsol[iCol] = 0.0; |
1397 | iCol = nextSlack[iCol]; |
1398 | } |
1399 | iCol = negSlack[i]; |
1400 | while (iCol >= 0) { |
1401 | colsol[iCol] = 0.0; |
1402 | iCol = nextSlack[iCol]; |
1403 | } |
1404 | { |
1405 | int iCol; |
1406 | iCol = posSlack[i]; |
1407 | while (iCol >= 0) { |
1408 | //printf("col %d el %g sol %g bounds %g %g cost %g\n", |
1409 | // iCol,element[columnStart[iCol]], |
1410 | // colsol[iCol],lower[iCol],upper[iCol],cost[iCol]); |
1411 | iCol = nextSlack[iCol]; |
1412 | } |
1413 | iCol = negSlack[i]; |
1414 | while (iCol >= 0) { |
1415 | //printf("col %d el %g sol %g bounds %g %g cost %g\n", |
1416 | // iCol,element[columnStart[iCol]], |
1417 | // colsol[iCol],lower[iCol],upper[iCol],cost[iCol]); |
1418 | iCol = nextSlack[iCol]; |
1419 | } |
1420 | } |
1421 | //printf("now what?\n"); |
1422 | int n = 0; |
1423 | bool basic = false; |
1424 | if (lo > 0.0) { |
1425 | // Add in positive |
1426 | iCol = posSlack[i]; |
1427 | while (iCol >= 0) { |
1428 | double value = element[columnStart[iCol]]; |
1429 | ratio[n] = cost[iCol] / value; |
1430 | which[n++] = iCol; |
1431 | iCol = nextSlack[iCol]; |
1432 | } |
1433 | CoinSort_2(ratio, ratio + n, which); |
1434 | for (int i = 0; i < n; i++) { |
1435 | iCol = which[i]; |
1436 | double value = element[columnStart[iCol]]; |
1437 | if (lo >= upper[iCol]*value) { |
1438 | value *= upper[iCol]; |
1439 | sum += value; |
1440 | lo -= value; |
1441 | colsol[iCol] = upper[iCol]; |
1442 | } else { |
1443 | value = lo / value; |
1444 | sum += lo; |
1445 | lo = 0.0; |
1446 | colsol[iCol] = value; |
1447 | model_->setColumnStatus(iCol, ClpSimplex::basic); |
1448 | basic = true; |
1449 | } |
1450 | if (lo < 1.0e-7) |
1451 | break; |
1452 | } |
1453 | } else if (up < 0.0) { |
1454 | // Use lo so coding is more similar |
1455 | lo = -up; |
1456 | // Add in negative |
1457 | iCol = negSlack[i]; |
1458 | while (iCol >= 0) { |
1459 | double value = -element[columnStart[iCol]]; |
1460 | ratio[n] = cost[iCol] / value; |
1461 | which[n++] = iCol; |
1462 | iCol = nextSlack[iCol]; |
1463 | } |
1464 | CoinSort_2(ratio, ratio + n, which); |
1465 | for (int i = 0; i < n; i++) { |
1466 | iCol = which[i]; |
1467 | double value = -element[columnStart[iCol]]; |
1468 | if (lo >= upper[iCol]*value) { |
1469 | value *= upper[iCol]; |
1470 | sum += value; |
1471 | lo -= value; |
1472 | colsol[iCol] = upper[iCol]; |
1473 | } else { |
1474 | value = lo / value; |
1475 | sum += lo; |
1476 | lo = 0.0; |
1477 | colsol[iCol] = value; |
1478 | model_->setColumnStatus(iCol, ClpSimplex::basic); |
1479 | basic = true; |
1480 | } |
1481 | if (lo < 1.0e-7) |
1482 | break; |
1483 | } |
1484 | } |
1485 | if (0) { |
1486 | double sum2 = 0.0; |
1487 | double x = 0.0; |
1488 | for (int k = 0; k < ncols; k++) { |
1489 | CoinBigIndex j; |
1490 | double value = colsol[k]; |
1491 | x += value * cost[k]; |
1492 | for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) { |
1493 | int irow = row[j]; |
1494 | if (irow == i) |
1495 | sum2 += element[j] * value; |
1496 | } |
1497 | } |
1498 | printf("after sum %g <= %g <= %g cost %.18g (sum = %g)\n" , |
1499 | rowlower[i], sum2, rowupper[i], x, sum); |
1500 | } |
1501 | rowsol[i] = sum; |
1502 | if (basic) { |
1503 | if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i])) |
1504 | model_->setRowStatus(i, ClpSimplex::atLowerBound); |
1505 | else |
1506 | model_->setRowStatus(i, ClpSimplex::atUpperBound); |
1507 | } |
1508 | } else { |
1509 | int n = 0; |
1510 | int iCol; |
1511 | iCol = posSlack[i]; |
1512 | while (iCol >= 0) { |
1513 | if (colsol[iCol] > lower[iCol] + 1.0e-8 && |
1514 | colsol[iCol] < upper[iCol] - 1.0e-8) { |
1515 | model_->setColumnStatus(iCol, ClpSimplex::basic); |
1516 | n++; |
1517 | } |
1518 | iCol = nextSlack[iCol]; |
1519 | } |
1520 | iCol = negSlack[i]; |
1521 | while (iCol >= 0) { |
1522 | if (colsol[iCol] > lower[iCol] + 1.0e-8 && |
1523 | colsol[iCol] < upper[iCol] - 1.0e-8) { |
1524 | model_->setColumnStatus(iCol, ClpSimplex::basic); |
1525 | n++; |
1526 | } |
1527 | iCol = nextSlack[iCol]; |
1528 | } |
1529 | if (n) { |
1530 | if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i])) |
1531 | model_->setRowStatus(i, ClpSimplex::atLowerBound); |
1532 | else |
1533 | model_->setRowStatus(i, ClpSimplex::atUpperBound); |
1534 | #ifdef CLP_INVESTIGATE |
1535 | if (n > 1) |
1536 | printf("%d basic on row %d!\n" , n, i); |
1537 | #endif |
1538 | } |
1539 | } |
1540 | } |
1541 | } |
1542 | delete [] ratio; |
1543 | delete [] which; |
1544 | #else |
1545 | for (i = 0; i < nrows; i++) { |
1546 | int n = 0; |
1547 | int iCol; |
1548 | iCol = posSlack[i]; |
1549 | while (iCol >= 0) { |
1550 | if (colsol[iCol] > lower[iCol] + 1.0e-8 && |
1551 | colsol[iCol] < upper[iCol] - 1.0e-8) { |
1552 | model_->setColumnStatus(iCol, ClpSimplex::basic); |
1553 | n++; |
1554 | } |
1555 | iCol = nextSlack[iCol]; |
1556 | } |
1557 | iCol = negSlack[i]; |
1558 | while (iCol >= 0) { |
1559 | if (colsol[iCol] > lower[iCol] + 1.0e-8 && |
1560 | colsol[iCol] < upper[iCol] - 1.0e-8) { |
1561 | model_->setColumnStatus(iCol, ClpSimplex::basic); |
1562 | n++; |
1563 | } |
1564 | iCol = nextSlack[iCol]; |
1565 | } |
1566 | if (n) { |
1567 | if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i])) |
1568 | model_->setRowStatus(i, ClpSimplex::atLowerBound); |
1569 | else |
1570 | model_->setRowStatus(i, ClpSimplex::atUpperBound); |
1571 | #ifdef CLP_INVESTIGATE |
1572 | if (n > 1) |
1573 | printf("%d basic on row %d!\n" , n, i); |
1574 | #endif |
1575 | } |
1576 | } |
1577 | #endif |
1578 | } |
1579 | double maxmin; |
1580 | if (model_->getObjSense() == -1.0) { |
1581 | maxmin = -1.0; |
1582 | } else { |
1583 | maxmin = 1.0; |
1584 | } |
1585 | bool justValuesPass = majorIterations_ > 1000000; |
1586 | if (slackStart >= 0) { |
1587 | for (i = 0; i < nrows; i++) { |
1588 | model_->setRowStatus(i, ClpSimplex::superBasic); |
1589 | } |
1590 | for (i = slackStart; i < slackEnd; i++) { |
1591 | model_->setColumnStatus(i, ClpSimplex::basic); |
1592 | } |
1593 | } else { |
1594 | /* still try and put singletons rather than artificials in basis */ |
1595 | int ninbas = 0; |
1596 | for (i = 0; i < nrows; i++) { |
1597 | model_->setRowStatus(i, ClpSimplex::basic); |
1598 | } |
1599 | for (i = 0; i < ncols; i++) { |
1600 | if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) { |
1601 | CoinBigIndex j = columnStart[i]; |
1602 | double value = element[j]; |
1603 | int irow = row[j]; |
1604 | double rlo = rowlower[irow]; |
1605 | double rup = rowupper[irow]; |
1606 | double clo = lower[i]; |
1607 | double cup = upper[i]; |
1608 | double csol = colsol[i]; |
1609 | /* adjust towards feasibility */ |
1610 | double move = 0.0; |
1611 | if (rowsol[irow] > rup) { |
1612 | move = (rup - rowsol[irow]) / value; |
1613 | if (value > 0.0) { |
1614 | /* reduce */ |
1615 | if (csol + move < clo) move = CoinMin(0.0, clo - csol); |
1616 | } else { |
1617 | /* increase */ |
1618 | if (csol + move > cup) move = CoinMax(0.0, cup - csol); |
1619 | } |
1620 | } else if (rowsol[irow] < rlo) { |
1621 | move = (rlo - rowsol[irow]) / value; |
1622 | if (value > 0.0) { |
1623 | /* increase */ |
1624 | if (csol + move > cup) move = CoinMax(0.0, cup - csol); |
1625 | } else { |
1626 | /* reduce */ |
1627 | if (csol + move < clo) move = CoinMin(0.0, clo - csol); |
1628 | } |
1629 | } else { |
1630 | /* move to improve objective */ |
1631 | if (cost[i]*maxmin > 0.0) { |
1632 | if (value > 0.0) { |
1633 | move = (rlo - rowsol[irow]) / value; |
1634 | /* reduce */ |
1635 | if (csol + move < clo) move = CoinMin(0.0, clo - csol); |
1636 | } else { |
1637 | move = (rup - rowsol[irow]) / value; |
1638 | /* increase */ |
1639 | if (csol + move > cup) move = CoinMax(0.0, cup - csol); |
1640 | } |
1641 | } else if (cost[i]*maxmin < 0.0) { |
1642 | if (value > 0.0) { |
1643 | move = (rup - rowsol[irow]) / value; |
1644 | /* increase */ |
1645 | if (csol + move > cup) move = CoinMax(0.0, cup - csol); |
1646 | } else { |
1647 | move = (rlo - rowsol[irow]) / value; |
1648 | /* reduce */ |
1649 | if (csol + move < clo) move = CoinMin(0.0, clo - csol); |
1650 | } |
1651 | } |
1652 | } |
1653 | rowsol[irow] += move * value; |
1654 | colsol[i] += move; |
1655 | /* put in basis if row was artificial */ |
1656 | if (rup - rlo < 1.0e-7 && model_->getRowStatus(irow) == ClpSimplex::basic) { |
1657 | model_->setRowStatus(irow, ClpSimplex::superBasic); |
1658 | model_->setColumnStatus(i, ClpSimplex::basic); |
1659 | ninbas++; |
1660 | } |
1661 | } |
1662 | } |
1663 | /*printf("%d in basis\n",ninbas);*/ |
1664 | } |
1665 | bool wantVector = false; |
1666 | if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) { |
1667 | // See if original wanted vector |
1668 | ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix()); |
1669 | wantVector = clpMatrixO->wantsSpecialColumnCopy(); |
1670 | } |
1671 | if (addAll < 3) { |
1672 | ClpPresolve pinfo; |
1673 | if (presolve) { |
1674 | if (allowInfeasible) { |
1675 | // fix up so will be feasible |
1676 | double * rhs = new double[nrows]; |
1677 | memset(rhs, 0, nrows * sizeof(double)); |
1678 | model_->clpMatrix()->times(1.0, colsol, rhs); |
1679 | double * rowupper = model_->rowUpper(); |
1680 | double * rowlower = model_->rowLower(); |
1681 | saveRowUpper = CoinCopyOfArray(rowupper, nrows); |
1682 | saveRowLower = CoinCopyOfArray(rowlower, nrows); |
1683 | double sum = 0.0; |
1684 | for (i = 0; i < nrows; i++) { |
1685 | if (rhs[i] > rowupper[i]) { |
1686 | sum += rhs[i] - rowupper[i]; |
1687 | rowupper[i] = rhs[i]; |
1688 | } |
1689 | if (rhs[i] < rowlower[i]) { |
1690 | sum += rowlower[i] - rhs[i]; |
1691 | rowlower[i] = rhs[i]; |
1692 | } |
1693 | } |
1694 | COIN_DETAIL_PRINT(printf("sum of infeasibilities %g\n" , sum)); |
1695 | delete [] rhs; |
1696 | } |
1697 | saveModel = model_; |
1698 | pinfo.setPresolveActions(pinfo.presolveActions() | 16384); |
1699 | model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5); |
1700 | } |
1701 | if (model_) { |
1702 | if (!wantVector) { |
1703 | //#define TWO_GOES |
1704 | #ifndef TWO_GOES |
1705 | model_->primal(justValuesPass ? 2 : 1); |
1706 | #else |
1707 | model_->primal(1 + 11); |
1708 | #endif |
1709 | } else { |
1710 | ClpMatrixBase * matrix = model_->clpMatrix(); |
1711 | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
1712 | assert (clpMatrix); |
1713 | clpMatrix->makeSpecialColumnCopy(); |
1714 | model_->primal(1); |
1715 | clpMatrix->releaseSpecialColumnCopy(); |
1716 | } |
1717 | if (presolve) { |
1718 | pinfo.postsolve(true); |
1719 | delete model_; |
1720 | model_ = saveModel; |
1721 | saveModel = NULL; |
1722 | } |
1723 | } else { |
1724 | // not feasible |
1725 | addAll = 1; |
1726 | presolve = 0; |
1727 | model_ = saveModel; |
1728 | saveModel = NULL; |
1729 | if (justValuesPass) |
1730 | model_->primal(2); |
1731 | } |
1732 | if (allowInfeasible) { |
1733 | CoinMemcpyN(saveRowUpper, nrows, model_->rowUpper()); |
1734 | CoinMemcpyN(saveRowLower, nrows, model_->rowLower()); |
1735 | delete [] saveRowUpper; |
1736 | delete [] saveRowLower; |
1737 | saveRowUpper = NULL; |
1738 | saveRowLower = NULL; |
1739 | } |
1740 | if (addAll < 2) { |
1741 | n = 0; |
1742 | if (!addAll ) { |
1743 | /* could do scans to get a good number */ |
1744 | iteration = 1; |
1745 | for (i = ordStart; i < ordEnd; i++) { |
1746 | if (whenUsed[i] >= iteration) { |
1747 | if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) { |
1748 | n++; |
1749 | upper[i] = saveUpper[i]; |
1750 | lower[i] = saveLower[i]; |
1751 | } |
1752 | } |
1753 | } |
1754 | } else { |
1755 | for (i = ordStart; i < ordEnd; i++) { |
1756 | if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) { |
1757 | n++; |
1758 | upper[i] = saveUpper[i]; |
1759 | lower[i] = saveLower[i]; |
1760 | } |
1761 | } |
1762 | delete [] saveUpper; |
1763 | delete [] saveLower; |
1764 | saveUpper = NULL; |
1765 | saveLower = NULL; |
1766 | } |
1767 | #ifdef COIN_DEVELOP |
1768 | printf("Time so far %g, %d now added from previous iterations\n" , |
1769 | CoinCpuTime() - startTime, n); |
1770 | #endif |
1771 | if (justValuesPass) |
1772 | return; |
1773 | if (addAll) |
1774 | presolve = 0; |
1775 | if (presolve) { |
1776 | saveModel = model_; |
1777 | model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5); |
1778 | } else { |
1779 | presolve = 0; |
1780 | } |
1781 | if (!wantVector) { |
1782 | model_->primal(1); |
1783 | } else { |
1784 | ClpMatrixBase * matrix = model_->clpMatrix(); |
1785 | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
1786 | assert (clpMatrix); |
1787 | clpMatrix->makeSpecialColumnCopy(); |
1788 | model_->primal(1); |
1789 | clpMatrix->releaseSpecialColumnCopy(); |
1790 | } |
1791 | if (presolve) { |
1792 | pinfo.postsolve(true); |
1793 | delete model_; |
1794 | model_ = saveModel; |
1795 | saveModel = NULL; |
1796 | } |
1797 | if (!addAll) { |
1798 | n = 0; |
1799 | for (i = ordStart; i < ordEnd; i++) { |
1800 | if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) { |
1801 | n++; |
1802 | upper[i] = saveUpper[i]; |
1803 | lower[i] = saveLower[i]; |
1804 | } |
1805 | } |
1806 | delete [] saveUpper; |
1807 | delete [] saveLower; |
1808 | saveUpper = NULL; |
1809 | saveLower = NULL; |
1810 | #ifdef COIN_DEVELOP |
1811 | printf("Time so far %g, %d now added from previous iterations\n" , |
1812 | CoinCpuTime() - startTime, n); |
1813 | #endif |
1814 | } |
1815 | if (presolve) { |
1816 | saveModel = model_; |
1817 | model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5); |
1818 | } else { |
1819 | presolve = 0; |
1820 | } |
1821 | if (!wantVector) { |
1822 | model_->primal(1); |
1823 | } else { |
1824 | ClpMatrixBase * matrix = model_->clpMatrix(); |
1825 | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
1826 | assert (clpMatrix); |
1827 | clpMatrix->makeSpecialColumnCopy(); |
1828 | model_->primal(1); |
1829 | clpMatrix->releaseSpecialColumnCopy(); |
1830 | } |
1831 | if (presolve) { |
1832 | pinfo.postsolve(true); |
1833 | delete model_; |
1834 | model_ = saveModel; |
1835 | saveModel = NULL; |
1836 | } |
1837 | } |
1838 | #ifdef COIN_DEVELOP |
1839 | printf("Total time in crossover %g\n" , CoinCpuTime() - startTime); |
1840 | #endif |
1841 | delete [] saveUpper; |
1842 | delete [] saveLower; |
1843 | } |
1844 | #ifdef FEB_TRY |
1845 | model_->setSpecialOptions(saveOptions); |
1846 | model_->setPerturbation(savePerturbation); |
1847 | #endif |
1848 | return ; |
1849 | } |
1850 | #endif |
1851 | /*****************************************************************************/ |
1852 | |
1853 | // Default contructor |
1854 | Idiot::Idiot() |
1855 | { |
1856 | model_ = NULL; |
1857 | maxBigIts_ = 3; |
1858 | maxIts_ = 5; |
1859 | logLevel_ = 1; |
1860 | logFreq_ = 100; |
1861 | maxIts2_ = 100; |
1862 | djTolerance_ = 1e-1; |
1863 | mu_ = 1e-4; |
1864 | drop_ = 5.0; |
1865 | exitDrop_ = -1.0e20; |
1866 | muFactor_ = 0.3333; |
1867 | stopMu_ = 1e-12; |
1868 | smallInfeas_ = 1e-1; |
1869 | reasonableInfeas_ = 1e2; |
1870 | muAtExit_ = 1.0e31; |
1871 | strategy_ = 8; |
1872 | lambdaIterations_ = 0; |
1873 | checkFrequency_ = 100; |
1874 | whenUsed_ = NULL; |
1875 | majorIterations_ = 30; |
1876 | exitFeasibility_ = -1.0; |
1877 | dropEnoughFeasibility_ = 0.02; |
1878 | dropEnoughWeighted_ = 0.01; |
1879 | // adjust |
1880 | double nrows = 10000.0; |
1881 | int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows))); |
1882 | baseIts = baseIts / 10; |
1883 | baseIts *= 10; |
1884 | maxIts2_ = 200 + baseIts + 5; |
1885 | maxIts2_ = 100; |
1886 | reasonableInfeas_ = static_cast<double> (nrows) * 0.05; |
1887 | lightWeight_ = 0; |
1888 | } |
1889 | // Constructor from model |
1890 | Idiot::Idiot(OsiSolverInterface &model) |
1891 | { |
1892 | model_ = & model; |
1893 | maxBigIts_ = 3; |
1894 | maxIts_ = 5; |
1895 | logLevel_ = 1; |
1896 | logFreq_ = 100; |
1897 | maxIts2_ = 100; |
1898 | djTolerance_ = 1e-1; |
1899 | mu_ = 1e-4; |
1900 | drop_ = 5.0; |
1901 | exitDrop_ = -1.0e20; |
1902 | muFactor_ = 0.3333; |
1903 | stopMu_ = 1e-12; |
1904 | smallInfeas_ = 1e-1; |
1905 | reasonableInfeas_ = 1e2; |
1906 | muAtExit_ = 1.0e31; |
1907 | strategy_ = 8; |
1908 | lambdaIterations_ = 0; |
1909 | checkFrequency_ = 100; |
1910 | whenUsed_ = NULL; |
1911 | majorIterations_ = 30; |
1912 | exitFeasibility_ = -1.0; |
1913 | dropEnoughFeasibility_ = 0.02; |
1914 | dropEnoughWeighted_ = 0.01; |
1915 | // adjust |
1916 | double nrows; |
1917 | if (model_) |
1918 | nrows = model_->getNumRows(); |
1919 | else |
1920 | nrows = 10000.0; |
1921 | int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows))); |
1922 | baseIts = baseIts / 10; |
1923 | baseIts *= 10; |
1924 | maxIts2_ = 200 + baseIts + 5; |
1925 | maxIts2_ = 100; |
1926 | reasonableInfeas_ = static_cast<double> (nrows) * 0.05; |
1927 | lightWeight_ = 0; |
1928 | } |
1929 | // Copy constructor. |
1930 | Idiot::Idiot(const Idiot &rhs) |
1931 | { |
1932 | model_ = rhs.model_; |
1933 | if (model_ && rhs.whenUsed_) { |
1934 | int numberColumns = model_->getNumCols(); |
1935 | whenUsed_ = new int [numberColumns]; |
1936 | CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_); |
1937 | } else { |
1938 | whenUsed_ = NULL; |
1939 | } |
1940 | djTolerance_ = rhs.djTolerance_; |
1941 | mu_ = rhs.mu_; |
1942 | drop_ = rhs.drop_; |
1943 | muFactor_ = rhs.muFactor_; |
1944 | stopMu_ = rhs.stopMu_; |
1945 | smallInfeas_ = rhs.smallInfeas_; |
1946 | reasonableInfeas_ = rhs.reasonableInfeas_; |
1947 | exitDrop_ = rhs.exitDrop_; |
1948 | muAtExit_ = rhs.muAtExit_; |
1949 | exitFeasibility_ = rhs.exitFeasibility_; |
1950 | dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_; |
1951 | dropEnoughWeighted_ = rhs.dropEnoughWeighted_; |
1952 | maxBigIts_ = rhs.maxBigIts_; |
1953 | maxIts_ = rhs.maxIts_; |
1954 | majorIterations_ = rhs.majorIterations_; |
1955 | logLevel_ = rhs.logLevel_; |
1956 | logFreq_ = rhs.logFreq_; |
1957 | checkFrequency_ = rhs.checkFrequency_; |
1958 | lambdaIterations_ = rhs.lambdaIterations_; |
1959 | maxIts2_ = rhs.maxIts2_; |
1960 | strategy_ = rhs.strategy_; |
1961 | lightWeight_ = rhs.lightWeight_; |
1962 | } |
1963 | // Assignment operator. This copies the data |
1964 | Idiot & |
1965 | Idiot::operator=(const Idiot & rhs) |
1966 | { |
1967 | if (this != &rhs) { |
1968 | delete [] whenUsed_; |
1969 | model_ = rhs.model_; |
1970 | if (model_ && rhs.whenUsed_) { |
1971 | int numberColumns = model_->getNumCols(); |
1972 | whenUsed_ = new int [numberColumns]; |
1973 | CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_); |
1974 | } else { |
1975 | whenUsed_ = NULL; |
1976 | } |
1977 | djTolerance_ = rhs.djTolerance_; |
1978 | mu_ = rhs.mu_; |
1979 | drop_ = rhs.drop_; |
1980 | muFactor_ = rhs.muFactor_; |
1981 | stopMu_ = rhs.stopMu_; |
1982 | smallInfeas_ = rhs.smallInfeas_; |
1983 | reasonableInfeas_ = rhs.reasonableInfeas_; |
1984 | exitDrop_ = rhs.exitDrop_; |
1985 | muAtExit_ = rhs.muAtExit_; |
1986 | exitFeasibility_ = rhs.exitFeasibility_; |
1987 | dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_; |
1988 | dropEnoughWeighted_ = rhs.dropEnoughWeighted_; |
1989 | maxBigIts_ = rhs.maxBigIts_; |
1990 | maxIts_ = rhs.maxIts_; |
1991 | majorIterations_ = rhs.majorIterations_; |
1992 | logLevel_ = rhs.logLevel_; |
1993 | logFreq_ = rhs.logFreq_; |
1994 | checkFrequency_ = rhs.checkFrequency_; |
1995 | lambdaIterations_ = rhs.lambdaIterations_; |
1996 | maxIts2_ = rhs.maxIts2_; |
1997 | strategy_ = rhs.strategy_; |
1998 | lightWeight_ = rhs.lightWeight_; |
1999 | } |
2000 | return *this; |
2001 | } |
2002 | Idiot::~Idiot() |
2003 | { |
2004 | delete [] whenUsed_; |
2005 | } |
2006 | |