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