source: stable/1.16/Clp/src/Idiot.cpp @ 2173

Last change on this file since 2173 was 2173, checked in by forrest, 4 years ago

out idiot assert

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 90.8 KB
Line 
1/* $Id: Idiot.cpp 2173 2015-09-30 10:24:02Z forrest $ */
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:
26For first two major iterations these are small.  Then:
27
28drop - exit a major iteration if drop over 5*checkFrequency < this is
29used as info->drop*(10.0+fabs(last weighted objective))
30
31exitDrop - exit idiot if feasible and drop < this is
32used as info->exitDrop*(10.0+fabs(last objective))
33
34djExit - exit a major iteration if largest dj (averaged over 5 checks)
35drops below this - used as info->djTolerance*(10.0+fabs(last weighted objective)
36
37djFlag - mostly skip variables with bad dj worse than this => 2*djExit
38
39djTol - 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
44int
45Idiot::dropping(IdiotResult result,
46                double tolerance,
47                double small,
48                int *nbad)
49{
50     if (result.infeas <= small) {
51          double value = CoinMax(fabs(result.objval), fabs(result.dropThis)) + 1.0;
52          if (result.dropThis > tolerance * value) {
53               *nbad = 0;
54               return 1;
55          } else {
56               (*nbad)++;
57               if (*nbad > 4) {
58                    return 0;
59               } else {
60                    return 1;
61               }
62          }
63     } else {
64          *nbad = 0;
65          return 1;
66     }
67}
68// Deals with whenUsed and slacks
69int
70Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd,
71                      double * colsol, const double * lower, const double * upper,
72                      const double * rowLower, const double * rowUpper,
73                      const double * cost, const double * element, double fixTolerance,
74                      double & objValue, double & infValue, double & maxInfeasibility)
75{
76     int n = 0;
77     if ((strategy_ & 16384) == 0) {
78          for (int i = ordinaryStart; i < ordinaryEnd; i++) {
79               if (colsol[i] > lower[i] + fixTolerance) {
80                    if (colsol[i] < upper[i] - fixTolerance) {
81                         n++;
82                    } else {
83                         colsol[i] = upper[i];
84                    }
85                    whenUsed_[i] = iteration;
86               } else {
87                    colsol[i] = lower[i];
88               }
89          }
90          return n;
91     } else {
92#ifdef COIN_DEVELOP
93          printf("entering inf %g, obj %g\n", infValue, objValue);
94#endif
95          int nrows = model_->getNumRows();
96          int ncols = model_->getNumCols();
97          int * posSlack = whenUsed_ + ncols;
98          int * negSlack = posSlack + nrows;
99          int * nextSlack = negSlack + nrows;
100          double * rowsol = reinterpret_cast<double *> (nextSlack + ncols);
101          memset(rowsol, 0, nrows * sizeof(double));
102#ifdef OSI_IDIOT
103          const CoinPackedMatrix * matrix = model_->getMatrixByCol();
104#else
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 &&
204                                   colsol[iCol] < upper[iCol] - fixTolerance) {
205                              whenUsed_[i] = iteration;
206                              n++;
207                         }
208                    }
209                    rowsol[i] = rowValue;
210               }
211               iCol = negSlack[i];
212               if (iCol >= 0) {
213                    // slide all slack down
214                    double rowValue = rowsol[i];
215                    CoinBigIndex j = columnStart[iCol];
216                    double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
217                    rowSave += (colsol[iCol] - lowerValue) * element[j];
218                    colsol[iCol] = lowerValue;
219                    while (nextSlack[iCol] >= 0) {
220                         iCol = nextSlack[iCol];
221                         j = columnStart[iCol];
222                         double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
223                         rowSave += (colsol[iCol] - lowerValue) * element[j];
224                         colsol[iCol] = lowerValue;
225                    }
226                    iCol = negSlack[i];
227                    while (rowValue > rowUpper[i] && iCol >= 0) {
228                         // want to increase
229                         double distance = -(rowUpper[i] - rowValue);
230                         double value = -element[columnStart[iCol]];
231                         double thisCost = cost[iCol];
232                         if (distance <= value*(upper[iCol] - lower[iCol])) {
233                              // can get there
234                              double movement = distance / value;
235                              objValue += movement * thisCost;
236                              rowValue = rowUpper[i];
237                              colsol[iCol] += movement;
238                         } else {
239                              // can't get there
240                              double movement = upper[iCol] - lower[iCol];
241                              objValue += movement * thisCost;
242                              rowValue -= movement * value;
243                              colsol[iCol] = upper[iCol];
244                              iCol = nextSlack[iCol];
245                         }
246                    }
247                    if (iCol >= 0) {
248                         // may want to carry on - because of cost?
249                         while (iCol >= 0 && cost[iCol] < 0 && rowValue > rowLower[i]) {
250                              // want to increase
251                              double distance = -(rowLower[i] - rowValue);
252                              double value = -element[columnStart[iCol]];
253                              double thisCost = cost[iCol];
254                              if (distance <= value*(upper[iCol] - colsol[iCol])) {
255                                   // can get there
256                                   double movement = distance / value;
257                                   objValue += movement * thisCost;
258                                   rowValue = rowLower[i];
259                                   colsol[iCol] += movement;
260                                   iCol = -1;
261                              } else {
262                                   // can't get there
263                                   double movement = upper[iCol] - colsol[iCol];
264                                   objValue += movement * thisCost;
265                                   rowValue -= movement * value;
266                                   colsol[iCol] = upper[iCol];
267                                   iCol = nextSlack[iCol];
268                              }
269                         }
270                         if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance &&
271                                   colsol[iCol] < upper[iCol] - fixTolerance) {
272                              whenUsed_[i] = iteration;
273                              n++;
274                         }
275                    }
276                    rowsol[i] = rowValue;
277               }
278               double infeasibility = CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]);
279               infValue += infeasibility;
280               maxInfeasibility=CoinMax(maxInfeasibility,infeasibility);
281               // just change
282               rowsol[i] -= rowSave;
283          }
284          return n;
285     }
286}
287
288/* returns -1 if none or start of costed slacks or -2 if
289   there are costed slacks but it is messy */
290static int countCostedSlacks(OsiSolverInterface * model)
291{
292#ifdef OSI_IDIOT
293     const CoinPackedMatrix * matrix = model->getMatrixByCol();
294#else
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) return -1;
309     while (1) {
310          for (i = 0; i < nrows; i++) {
311               int j = i + slackStart;
312               CoinBigIndex k = columnStart[j];
313               if (columnLength[j] == 1) {
314                    if (row[k] != i || element[k] != 1.0) {
315                         nSlacks = 0;
316                         break;
317                    }
318               } else {
319                    nSlacks = 0;
320                    break;
321               }
322               if (rowupper[i] <= 0.0) {
323                    nSlacks = 0;
324                    break;
325               }
326          }
327          if (nSlacks || !slackStart) break;
328          slackStart = 0;
329     }
330     if (!nSlacks) slackStart = -1;
331     return slackStart;
332}
333void
334Idiot::crash(int numberPass, CoinMessageHandler * handler,
335             const CoinMessages *messages, bool doCrossover)
336{
337     // lightweight options
338     int numberColumns = model_->getNumCols();
339     const double * objective = model_->getObjCoefficients();
340     int nnzero = 0;
341     double sum = 0.0;
342     int i;
343     for (i = 0; i < numberColumns; i++) {
344          if (objective[i]) {
345               sum += fabs(objective[i]);
346               nnzero++;
347          }
348     }
349     sum /= static_cast<double> (nnzero + 1);
350     if (maxIts_ == 5)
351          maxIts_ = 2;
352     if (numberPass <= 0)
353          majorIterations_ = static_cast<int>(2 + log10(static_cast<double>(numberColumns + 1)));
354     else
355          majorIterations_ = numberPass;
356     // If mu not changed then compute
357     if (mu_ == 1e-4)
358          mu_ = CoinMax(1.0e-3, sum * 1.0e-5);
359     if (maxIts2_ == 100) {
360          if (!lightWeight_) {
361               maxIts2_ = 105;
362          } else if (lightWeight_ == 1) {
363               mu_ *= 1000.0;
364               maxIts2_ = 23;
365          } else if (lightWeight_ == 2) {
366               maxIts2_ = 11;
367          } else {
368               maxIts2_ = 23;
369          }
370     }
371     //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
372     solve2(handler, messages);
373#ifndef OSI_IDIOT
374     if (doCrossover) {
375          double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows());
376          if ((averageInfeas < 0.01 && (strategy_ & 512) != 0) || (strategy_ & 8192) != 0)
377               crossOver(16 + 1);
378          else
379               crossOver(majorIterations_ < 1000000 ? 3 : 2);
380     }
381#endif
382}
383
384void
385Idiot::solve()
386{
387     CoinMessages dummy;
388     solve2(NULL, &dummy);
389}
390void
391Idiot::solve2(CoinMessageHandler * handler, const CoinMessages * messages)
392{
393     int strategy = 0;
394     double d2;
395     int i, n;
396     int allOnes = 1;
397     int iteration = 0;
398     int iterationTotal = 0;
399     int nTry = 0; /* number of tries at same weight */
400     double fixTolerance = IDIOT_FIX_TOLERANCE;
401     int maxBigIts = maxBigIts_;
402     int maxIts = maxIts_;
403     int logLevel = logLevel_;
404     int saveMajorIterations = majorIterations_;
405     majorIterations_ = majorIterations_ % 1000000;
406     if (handler) {
407          if (handler->logLevel() > 0 && handler->logLevel() < 3)
408               logLevel = 1;
409          else if (!handler->logLevel())
410               logLevel = 0;
411          else
412               logLevel = 7;
413     }
414     double djExit = djTolerance_;
415     double djFlag = 1.0 + 100.0 * djExit;
416     double djTol = 0.00001;
417     double mu = mu_;
418     double drop = drop_;
419     int maxIts2 = maxIts2_;
420     double factor = muFactor_;
421     double smallInfeas = smallInfeas_;
422     double reasonableInfeas = reasonableInfeas_;
423     double stopMu = stopMu_;
424     double maxmin, offset;
425     double lastWeighted = 1.0e50;
426     double exitDrop = exitDrop_;
427     double fakeSmall = smallInfeas;
428     double firstInfeas;
429     int badIts = 0;
430     int slackStart, ordStart, ordEnd;
431     int checkIteration = 0;
432     int lambdaIteration = 0;
433     int belowReasonable = 0; /* set if ever gone below reasonable infeas */
434     double bestWeighted = 1.0e60;
435     double bestFeasible = 1.0e60; /* best solution while feasible */
436     IdiotResult result, lastResult;
437     int saveStrategy = strategy_;
438     const int strategies[] = {0, 2, 128};
439     double saveLambdaScale = 0.0;
440     if ((saveStrategy & 128) != 0) {
441          fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
442     }
443#ifdef OSI_IDIOT
444     const CoinPackedMatrix * matrix = model_->getMatrixByCol();
445#else
446     ClpMatrixBase * matrix = model_->clpMatrix();
447#endif
448     const int * row = matrix->getIndices();
449     const CoinBigIndex * columnStart = matrix->getVectorStarts();
450     const int * columnLength = matrix->getVectorLengths();
451     const double * element = matrix->getElements();
452     int nrows = model_->getNumRows();
453     int ncols = model_->getNumCols();
454     double * rowsol, * colsol;
455     double * pi, * dj;
456#ifndef OSI_IDIOT
457     double * cost = model_->objective();
458     double * lower = model_->columnLower();
459     double * upper = model_->columnUpper();
460#else
461     double * cost = new double [ncols];
462     CoinMemcpyN( model_->getObjCoefficients(), ncols, cost);
463     const double * lower = model_->getColLower();
464     const double * upper = model_->getColUpper();
465#endif
466     const double *elemXX;
467     double * saveSol;
468     double * rowupper = new double[nrows]; // not const as modified
469     CoinMemcpyN(model_->getRowUpper(), nrows, rowupper);
470     double * rowlower = new double[nrows]; // not const as modified
471     CoinMemcpyN(model_->getRowLower(), nrows, rowlower);
472     CoinThreadRandom * randomNumberGenerator = model_->randomNumberGenerator();
473     int * whenUsed;
474     double * lambda;
475     saveSol = new double[ncols];
476     lambda = new double [nrows];
477     rowsol = new double[nrows];
478     colsol = new double [ncols];
479     CoinMemcpyN(model_->getColSolution(), ncols, colsol);
480     pi = new double[nrows];
481     dj = new double[ncols];
482#ifndef OSI_IDIOT
483     bool fixAfterSome = false; //(model_->specialOptions()&8388608)!=0;
484     int exitAfter = 50; //(model_->specialOptions()&8388608)!=0 ? 50 : 1000000;
485     {
486        int numberColumns = model_->numberColumns();
487        for (int i=0;i<numberColumns;i++) {
488          if (upper[i]==lower[i])
489            model_->setColumnStatus(i,ClpSimplex::isFixed);
490        }
491     }
492#endif
493     delete [] whenUsed_;
494     bool oddSlacks = false;
495     // See if any costed slacks
496     int numberSlacks = 0;
497     for (i = 0; i < ncols; i++) {
498          if (columnLength[i] == 1)
499               numberSlacks++;
500     }
501     if (!numberSlacks || (strategy_&524288)!=0) {
502          whenUsed_ = new int[ncols];
503     } else {
504#ifdef COIN_DEVELOP
505          printf("%d slacks\n", numberSlacks);
506#endif
507          oddSlacks = true;
508          int extra = static_cast<int> (nrows * sizeof(double) / sizeof(int));
509          whenUsed_ = new int[2*ncols+2*nrows+extra];
510          int * posSlack = whenUsed_ + ncols;
511          int * negSlack = posSlack + nrows;
512          int * nextSlack = negSlack + nrows;
513          for (i = 0; i < nrows; i++) {
514               posSlack[i] = -1;
515               negSlack[i] = -1;
516          }
517          for (i = 0; i < ncols; i++)
518               nextSlack[i] = -1;
519          for (i = 0; i < ncols; i++) {
520               if (columnLength[i] == 1) {
521                    CoinBigIndex j = columnStart[i];
522                    int iRow = row[j];
523                    if (element[j] > 0.0) {
524                         if (posSlack[iRow] == -1) {
525                              posSlack[iRow] = i;
526                         } else {
527                              int iCol = posSlack[iRow];
528                              while (nextSlack[iCol] >= 0)
529                                   iCol = nextSlack[iCol];
530                              nextSlack[iCol] = i;
531                         }
532                    } else {
533                         if (negSlack[iRow] == -1) {
534                              negSlack[iRow] = i;
535                         } else {
536                              int iCol = negSlack[iRow];
537                              while (nextSlack[iCol] >= 0)
538                                   iCol = nextSlack[iCol];
539                              nextSlack[iCol] = i;
540                         }
541                    }
542               }
543          }
544          // now sort
545          for (i = 0; i < nrows; i++) {
546               int iCol;
547               iCol = posSlack[i];
548               if (iCol >= 0) {
549                    CoinBigIndex j = columnStart[iCol];
550#ifndef NDEBUG
551                    int iRow = row[j];
552#endif
553                    assert (element[j] > 0.0);
554                    assert (iRow == i);
555                    dj[0] = cost[iCol] / element[j];
556                    whenUsed_[0] = iCol;
557                    int n = 1;
558                    while (nextSlack[iCol] >= 0) {
559                         iCol = nextSlack[iCol];
560                         CoinBigIndex j = columnStart[iCol];
561#ifndef NDEBUG
562                         int iRow = row[j];
563#endif
564                         assert (element[j] > 0.0);
565                         assert (iRow == i);
566                         dj[n] = cost[iCol] / element[j];
567                         whenUsed_[n++] = iCol;
568                    }
569                    for (j = 0; j < n; j++) {
570                         int jCol = whenUsed_[j];
571                         nextSlack[jCol] = -2;
572                    }
573                    CoinSort_2(dj, dj + n, whenUsed_);
574                    // put back
575                    iCol = whenUsed_[0];
576                    posSlack[i] = iCol;
577                    for (j = 1; j < n; j++) {
578                         int jCol = whenUsed_[j];
579                         nextSlack[iCol] = jCol;
580                         iCol = jCol;
581                    }
582               }
583               iCol = negSlack[i];
584               if (iCol >= 0) {
585                    CoinBigIndex j = columnStart[iCol];
586#ifndef NDEBUG
587                    int iRow = row[j];
588#endif
589                    assert (element[j] < 0.0);
590                    assert (iRow == i);
591                    dj[0] = -cost[iCol] / element[j];
592                    whenUsed_[0] = iCol;
593                    int n = 1;
594                    while (nextSlack[iCol] >= 0) {
595                         iCol = nextSlack[iCol];
596                         CoinBigIndex j = columnStart[iCol];
597#ifndef NDEBUG
598                         int iRow = row[j];
599#endif
600                         assert (element[j] < 0.0);
601                         assert (iRow == i);
602                         dj[n] = -cost[iCol] / element[j];
603                         whenUsed_[n++] = iCol;
604                    }
605                    for (j = 0; j < n; j++) {
606                         int jCol = whenUsed_[j];
607                         nextSlack[jCol] = -2;
608                    }
609                    CoinSort_2(dj, dj + n, whenUsed_);
610                    // put back
611                    iCol = whenUsed_[0];
612                    negSlack[i] = iCol;
613                    for (j = 1; j < n; j++) {
614                         int jCol = whenUsed_[j];
615                         nextSlack[iCol] = jCol;
616                         iCol = jCol;
617                    }
618               }
619          }
620     }
621     whenUsed = whenUsed_;
622     if (model_->getObjSense() == -1.0) {
623          maxmin = -1.0;
624     } else {
625          maxmin = 1.0;
626     }
627     model_->getDblParam(OsiObjOffset, offset);
628     if (!maxIts2) maxIts2 = maxIts;
629     strategy = strategy_;
630     strategy &= 3;
631     memset(lambda, 0, nrows * sizeof(double));
632     slackStart = countCostedSlacks(model_);
633     // redo in case odd matrix
634     row = matrix->getIndices();
635     columnStart = matrix->getVectorStarts();
636     columnLength = matrix->getVectorLengths();
637     element = matrix->getElements();
638     if (slackStart >= 0) {
639       COIN_DETAIL_PRINT(printf("This model has costed slacks\n"));
640          if (slackStart) {
641               ordStart = 0;
642               ordEnd = slackStart;
643          } else {
644               ordStart = nrows;
645               ordEnd = ncols;
646          }
647     } else {
648          ordStart = 0;
649          ordEnd = ncols;
650     }
651     if (offset && logLevel > 2) {
652          printf("** Objective offset is %g\n", offset);
653     }
654     /* compute reasonable solution cost */
655     for (i = 0; i < nrows; i++) {
656          rowsol[i] = 1.0e31;
657     }
658     for (i = 0; i < ncols; i++) {
659          CoinBigIndex j;
660          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
661               if (element[j] != 1.0) {
662                    allOnes = 0;
663                    break;
664               }
665          }
666     }
667     if (allOnes) {
668          elemXX = NULL;
669     } else {
670          elemXX = element;
671     }
672     // Do scaling if wanted
673     bool scaled = false;
674#ifndef OSI_IDIOT
675     if ((strategy_ & 32) != 0 && !allOnes) {
676          if (model_->scalingFlag() > 0)
677               scaled = model_->clpMatrix()->scale(model_) == 0;
678          if (scaled) {
679#define IDIOT_SCALE 2
680#ifndef IDIOT_SCALE
681               const double * rowScale = model_->rowScale();
682#else
683               double * rowScale = model_->mutableRowScale();
684#endif
685               const double * columnScale = model_->columnScale();
686               double * oldLower = lower;
687               double * oldUpper = upper;
688               double * oldCost = cost;
689               lower = new double[ncols];
690               upper = new double[ncols];
691               cost = new double[ncols];
692               CoinMemcpyN(oldLower, ncols, lower);
693               CoinMemcpyN(oldUpper, ncols, upper);
694               CoinMemcpyN(oldCost, ncols, cost);
695               int icol, irow;
696#if IDIOT_SCALE<0
697               for (irow = 0; irow < nrows; irow++) {
698                 rowlower[irow]=1.0e100;
699                 rowupper[irow]=1.0e-100;
700               }
701#endif
702               for (icol = 0; icol < ncols; icol++) {
703                    double multiplier = 1.0 / columnScale[icol];
704                    if (lower[icol] > -1.0e50)
705                         lower[icol] *= multiplier;
706                    if (upper[icol] < 1.0e50)
707                         upper[icol] *= multiplier;
708                    colsol[icol] *= multiplier;
709                    cost[icol] *= columnScale[icol];
710#if IDIOT_SCALE<0
711                    CoinBigIndex j;
712                    double scale = columnScale[i];
713                    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
714                         int jrow = row[j];
715                         double scaledValue = fabs(scale*element[j]);
716                         rowlower[jrow]=CoinMin(rowlower[jrow],scaledValue);
717                         rowupper[jrow]=CoinMax(rowupper[jrow],scaledValue);
718                    }
719#endif
720               }
721#ifdef IDIOT_SCALE
722#if IDIOT_SCALE>1||IDIOT_SCALE<-1
723               const double * rowLower=model_->rowLower();
724               const double * rowUpper=model_->rowUpper();
725#endif
726               for (irow = 0; irow < nrows; irow++) {
727#if IDIOT_SCALE<0
728                 double multiplier = 1.0/sqrt(rowlower[irow]*rowupper[irow]);
729#else
730                 double multiplier = rowScale[irow];
731#endif
732#if IDIOT_SCALE>1||IDIOT_SCALE<-1
733#define EQUALITY_MULTIPLIER 2
734                 if (rowLower[irow]==rowUpper[irow])
735                   multiplier *= EQUALITY_MULTIPLIER;
736#if IDIOT_SCALE>2||IDIOT_SCALE<-2
737                 if (rowLower[irow]==rowUpper[irow]&&!rowlower[irow])
738                   multiplier *= EQUALITY_MULTIPLIER;
739#endif
740#endif
741                 rowScale[irow]=multiplier;
742               }
743               CoinMemcpyN(model_->rowUpper(), nrows, rowupper);
744#endif
745               CoinMemcpyN(model_->rowLower(), nrows, rowlower);
746               for (irow = 0; irow < nrows; irow++) {
747                    double multiplier = rowScale[irow];
748                    if (rowlower[irow] > -1.0e50)
749                         rowlower[irow] *= multiplier;
750                    if (rowupper[irow] < 1.0e50)
751                         rowupper[irow] *= multiplier;
752                    rowsol[irow] *= multiplier;
753               }
754               int length = columnStart[ncols-1] + columnLength[ncols-1];
755               double * elemYY = new double[length];
756               for (i = 0; i < ncols; i++) {
757                    CoinBigIndex j;
758                    double scale = columnScale[i];
759                    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
760                         int irow = row[j];
761                         elemYY[j] = element[j] * scale * rowScale[irow];
762                    }
763               }
764               elemXX = elemYY;
765          }
766     }
767#endif
768     if ((strategy_&131072)!=0) {
769       // going to mess with cost
770       if (cost==model_->objective())
771         cost=CoinCopyOfArray(cost,ncols);
772     }
773     for (i = 0; i < ncols; i++) {
774          CoinBigIndex j;
775          double dd = columnLength[i];
776          dd = cost[i] / dd;
777          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
778               int irow = row[j];
779               if (dd < rowsol[irow]) {
780                    rowsol[irow] = dd;
781               }
782          }
783     }
784     d2 = 0.0;
785     for (i = 0; i < nrows; i++) {
786          d2 += rowsol[i];
787     }
788     d2 *= 2.0; /* for luck */
789
790     d2 = d2 / static_cast<double> (4 * nrows + 8000);
791     d2 *= 0.5; /* halve with more flexible method */
792     if (d2 < 5.0) d2 = 5.0;
793     if (djExit == 0.0) {
794          djExit = d2;
795     }
796     if ((saveStrategy & 4) != 0) {
797          /* go to relative tolerances - first small */
798          djExit = 1.0e-10;
799          djFlag = 1.0e-5;
800          drop = 1.0e-10;
801     }
802     memset(whenUsed, 0, ncols * sizeof(int));
803     strategy = strategies[strategy];
804     if ((saveStrategy & 8) != 0) strategy |= 64; /* don't allow large theta */
805     CoinMemcpyN(colsol, ncols, saveSol);
806
807     lastResult = IdiSolve(nrows, ncols, rowsol , colsol, pi,
808                           dj, cost, rowlower, rowupper,
809                           lower, upper, elemXX, row, columnStart, columnLength, lambda,
810                           0, mu, drop,
811                           maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
812     // update whenUsed_
813     double maxInfeasibility=COIN_DBL_MAX;
814     n = cleanIteration(iteration, ordStart, ordEnd,
815                        colsol,  lower,  upper,
816                        rowlower, rowupper,
817                        cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas, maxInfeasibility);
818     if ((strategy_ & 16384) != 0) {
819          int * posSlack = whenUsed_ + ncols;
820          int * negSlack = posSlack + nrows;
821          int * nextSlack = negSlack + nrows;
822          double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
823          for (i = 0; i < nrows; i++)
824               rowsol[i] += rowsol2[i];
825     }
826     if ((logLevel_ & 1) != 0) {
827#ifndef OSI_IDIOT
828          if (!handler) {
829#endif
830               printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
831                      iteration, lastResult.infeas, lastResult.objval, mu, lastResult.iteration, n);
832#ifndef OSI_IDIOT
833          } else {
834               handler->message(CLP_IDIOT_ITERATION, *messages)
835                         << iteration << lastResult.infeas << lastResult.objval << mu << lastResult.iteration << n
836                         << CoinMessageEol;
837          }
838#endif
839     }
840     int numberBaseTrys = 0; // for first time
841     int numberAway = -1;
842     iterationTotal = lastResult.iteration;
843     firstInfeas = lastResult.infeas;
844     if ((strategy_ & 1024) != 0) reasonableInfeas = 0.5 * firstInfeas;
845     if (lastResult.infeas < reasonableInfeas) lastResult.infeas = reasonableInfeas;
846     double keepinfeas = 1.0e31;
847     double lastInfeas = 1.0e31;
848     double bestInfeas = 1.0e31;
849     while ((mu > stopMu && lastResult.infeas > smallInfeas) ||
850               (lastResult.infeas <= smallInfeas &&
851                dropping(lastResult, exitDrop, smallInfeas, &badIts)) ||
852               checkIteration < 2 || lambdaIteration < lambdaIterations_) {
853          if (lastResult.infeas <= exitFeasibility_)
854               break;
855          iteration++;
856          //if (iteration>=exitAfter)
857          //break;
858          checkIteration++;
859          if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) {
860               bestFeasible = lastResult.objval;
861          }
862          if (lastResult.infeas + mu * lastResult.objval < bestWeighted) {
863               bestWeighted = lastResult.objval + mu * lastResult.objval;
864          }
865          if ((saveStrategy & 4096)) strategy |= 256;
866          if ((saveStrategy & 4) != 0 && iteration > 2) {
867               /* go to relative tolerances */
868               double weighted = 10.0 + fabs(lastWeighted);
869               djExit = djTolerance_ * weighted;
870               djFlag = 2.0 * djExit;
871               drop = drop_ * weighted;
872               djTol = 0.01 * djExit;
873          }
874          CoinMemcpyN(colsol, ncols, saveSol);
875          result = IdiSolve(nrows, ncols, rowsol , colsol, pi, dj,
876                            cost, rowlower, rowupper,
877                            lower, upper, elemXX, row, columnStart, columnLength, lambda,
878                            maxIts, mu, drop,
879                            maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
880          n = cleanIteration(iteration, ordStart, ordEnd,
881                             colsol,  lower,  upper,
882                             rowlower, rowupper,
883                             cost, elemXX, fixTolerance, result.objval, result.infeas, maxInfeasibility);
884          if ((strategy_ & 16384) != 0) {
885               int * posSlack = whenUsed_ + ncols;
886               int * negSlack = posSlack + nrows;
887               int * nextSlack = negSlack + nrows;
888               double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
889               for (i = 0; i < nrows; i++)
890                    rowsol[i] += rowsol2[i];
891          } else {
892            maxInfeasibility=0.0;
893            for (i = 0; i < nrows; i++) {
894              double infeasibility = CoinMax(CoinMax(0.0, rowlower[i] - rowsol[i]), rowsol[i] - rowupper[i]);
895              maxInfeasibility=CoinMax(maxInfeasibility,infeasibility);
896            }
897          }
898          if ((logLevel_ & 1) != 0) {
899#ifndef OSI_IDIOT
900               if (!handler) {
901#endif
902                    printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
903                           iteration, result.infeas, result.objval, mu, result.iteration, n);
904#ifndef OSI_IDIOT
905               } else {
906                    handler->message(CLP_IDIOT_ITERATION, *messages)
907                              << iteration << result.infeas << result.objval << mu << result.iteration << n
908                              << CoinMessageEol;
909               }
910#endif
911          }
912#ifndef OSI_IDIOT
913          if (fixAfterSome) {
914            if (result.infeas<0.01*nrows&&iteration>10&&(3*n>2*nrows||4*n>2*ncols)) {
915                // flag
916                int numberColumns = model_->numberColumns();
917                printf("Flagging satisfied\n");
918                fixAfterSome=false;
919                for (int i=0;i<numberColumns;i++) {
920                  if (colsol[i]>upper[i]-1.0e-7||
921                      colsol[i]<lower[i]+1.0e-7)
922                    model_->setColumnStatus(i,ClpSimplex::isFixed);
923                }
924              }
925              }
926#endif
927          if (iteration > exitAfter ) {
928            if((result.infeas < 1.0e-4 && majorIterations_<200 && n == numberAway)||result.infeas<1.0e-8||maxInfeasibility<1.0e-8) {
929#ifdef CLP_INVESTIGATE
930               printf("infeas small %g\n", result.infeas);
931#endif
932               if ((strategy_&131072)==0) {
933                 break; // not much happening
934               } else {
935                 int n=0;
936                 for (int i=0;i<ncols;i++) {
937                   if (cost[i])
938                     n++;
939                 }
940                 if (n*10<ncols) {
941                   // fix ones with costs
942                   exitAfter=100;
943                   for (int i=0;i<ncols;i++) {
944                     if (cost[i]||colsol[i]<lower[i]+1.0e-9||
945                         colsol[i]>upper[i]-1.0e-9) {
946                       cost[i]=0.0;
947                       model_->setColumnStatus(i,ClpSimplex::isFixed);
948                     } else {
949                       cost[i]=0.0;
950                       double gap = upper[i]-lower[i];
951                       if (colsol[i]>upper[i]-0.25*gap) {
952                         cost[i]=gap/(colsol[i]-upper[i]);
953                       } else if (colsol[i]<lower[i]+0.25*gap) {
954                         cost[i]=gap/(colsol[i]-lower[i]);
955                       }
956                     }
957                   }
958                 }
959               }
960            }
961          }
962          if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) {
963               if (lastInfeas != bestInfeas && CoinMin(result.infeas, lastInfeas) > 0.95 * bestInfeas)
964                    majorIterations_ = CoinMin(majorIterations_, iteration); // not getting feasible
965          }
966          lastInfeas = result.infeas;
967          numberAway = n;
968          keepinfeas = result.infeas;
969          lastWeighted = result.weighted;
970          iterationTotal += result.iteration;
971          if (iteration == 1) {
972               if ((strategy_ & 1024) != 0 && mu < 1.0e-10)
973                    result.infeas = firstInfeas * 0.8;
974               if (majorIterations_ >= 50 || dropEnoughFeasibility_ <= 0.0)
975                    result.infeas *= 0.8;
976               if (result.infeas > firstInfeas * 0.9
977                         && result.infeas > reasonableInfeas) {
978                    iteration--;
979                    if (majorIterations_ < 50)
980                         mu *= 1.0e-1;
981                    else
982                         mu *= 0.7;
983                    bestFeasible = 1.0e31;
984                    bestWeighted = 1.0e60;
985                    numberBaseTrys++;
986                    if (mu < 1.0e-30 || (numberBaseTrys > 10 && lightWeight_)) {
987                         // back to all slack basis
988                         lightWeight_ = 2;
989                         break;
990                    }
991                    CoinMemcpyN(saveSol, ncols, colsol);
992               } else {
993                    maxIts = maxIts2;
994                    checkIteration = 0;
995                    if ((strategy_ & 1024) != 0) mu *= 1.0e-1;
996               }
997          } else {
998          }
999          bestInfeas = CoinMin(bestInfeas, result.infeas);
1000          if (majorIterations_>100&&majorIterations_<200) {
1001            if (iteration==majorIterations_-100) {
1002              // redo
1003              double muX=mu*10.0;
1004              bestInfeas=1.0e3;
1005              mu=muX;
1006              nTry=0;
1007            }
1008          }
1009          if (iteration) {
1010               /* this code is in to force it to terminate sometime */
1011               double changeMu = factor;
1012               if ((saveStrategy & 64) != 0) {
1013                    keepinfeas = 0.0; /* switch off ranga's increase */
1014                    fakeSmall = smallInfeas;
1015               } else {
1016                    fakeSmall = -1.0;
1017               }
1018               saveLambdaScale = 0.0;
1019               if (result.infeas > reasonableInfeas ||
1020                         (nTry + 1 == maxBigIts && result.infeas > fakeSmall)) {
1021                    if (result.infeas > lastResult.infeas*(1.0 - dropEnoughFeasibility_) ||
1022                              nTry + 1 == maxBigIts ||
1023                              (result.infeas > lastResult.infeas * 0.9
1024                               && result.weighted > lastResult.weighted
1025                               - 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) drop = drop_ / 50.0;
1055                    if (lambdaIteration > 4 &&
1056                              (((lambdaIteration % 10) == 0 && smallInfeas < keepinfeas) ||
1057                               ((lambdaIteration % 5) == 0 && 1.5 * smallInfeas < keepinfeas))) {
1058                         //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
1059                         smallInfeas *= 1.5;
1060                    }
1061                    if ((saveStrategy & 2048) == 0) {
1062                         for (i = 0; i < nrows; i++) {
1063                              if (lambda[i]) 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
1297void
1298Idiot::crossOver(int mode)
1299{
1300     if (lightWeight_ == 2) {
1301          // total failure
1302          model_->allSlackBasis();
1303          return;
1304     }
1305     double fixTolerance = IDIOT_FIX_TOLERANCE;
1306#ifdef COIN_DEVELOP
1307     double startTime = CoinCpuTime();
1308#endif
1309     ClpSimplex * saveModel = NULL;
1310     ClpMatrixBase * matrix = model_->clpMatrix();
1311     const int * row = matrix->getIndices();
1312     const CoinBigIndex * columnStart = matrix->getVectorStarts();
1313     const int * columnLength = matrix->getVectorLengths();
1314     const double * element = matrix->getElements();
1315     const double * rowupper = model_->getRowUpper();
1316     model_->eventHandler()->event(ClpEventHandler::startOfCrossover);
1317     int nrows = model_->getNumRows();
1318     int ncols = model_->getNumCols();
1319     double * rowsol, * colsol;
1320     // different for Osi
1321     double * lower = model_->columnLower();
1322     double * upper = model_->columnUpper();
1323     const double * rowlower = model_->getRowLower();
1324     int * whenUsed = whenUsed_;
1325     rowsol = model_->primalRowSolution();
1326     colsol = model_->primalColumnSolution();;
1327     double * cost = model_->objective();
1328
1329     int slackEnd, ordStart, ordEnd;
1330     int slackStart = countCostedSlacks(model_);
1331
1332     int addAll = mode & 7;
1333     int presolve = 0;
1334
1335     double djTolerance = djTolerance_;
1336     if (djTolerance > 0.0 && djTolerance < 1.0)
1337          djTolerance = 1.0;
1338     int iteration;
1339     int i, n = 0;
1340     double ratio = 1.0;
1341     double objValue = 0.0;
1342     if ((strategy_ & 128) != 0) {
1343          fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
1344     }
1345     if ((mode & 16) != 0 && addAll < 3) 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 ||
1502                              rowupper[i] - rowsol[i] < 1.0e-7) &&
1503                              nPlus + nMinus < 2)
1504                         possible = false;
1505                    if (possible) {
1506                         // Amount contributed by other varaibles
1507                         sum = rowsol[i] - sum;
1508                         double lo = rowlower[i];
1509                         if (lo > -1.0e20)
1510                              lo -= sum;
1511                         double up = rowupper[i];
1512                         if (up < 1.0e20)
1513                              up -= sum;
1514                         //printf("row bounds %g %g\n",lo,up);
1515                         if (0) {
1516                              double sum = 0.0;
1517                              double x = 0.0;
1518                              for (int k = 0; k < ncols; k++) {
1519                                   CoinBigIndex j;
1520                                   double value = colsol[k];
1521                                   x += value * cost[k];
1522                                   for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
1523                                        int irow = row[j];
1524                                        if (irow == i)
1525                                             sum += element[j] * value;
1526                                   }
1527                              }
1528                              printf("Before sum %g <= %g <= %g cost %.18g\n",
1529                                     rowlower[i], sum, rowupper[i], x);
1530                         }
1531                         // set all to zero
1532                         iCol = posSlack[i];
1533                         while (iCol >= 0) {
1534                              colsol[iCol] = 0.0;
1535                              iCol = nextSlack[iCol];
1536                         }
1537                         iCol = negSlack[i];
1538                         while (iCol >= 0) {
1539                              colsol[iCol] = 0.0;
1540                              iCol = nextSlack[iCol];
1541                         }
1542                         {
1543                              int iCol;
1544                              iCol = posSlack[i];
1545                              while (iCol >= 0) {
1546                                   //printf("col %d el %g sol %g bounds %g %g cost %g\n",
1547                                   //     iCol,element[columnStart[iCol]],
1548                                   //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
1549                                   iCol = nextSlack[iCol];
1550                              }
1551                              iCol = negSlack[i];
1552                              while (iCol >= 0) {
1553                                   //printf("col %d el %g sol %g bounds %g %g cost %g\n",
1554                                   //     iCol,element[columnStart[iCol]],
1555                                   //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
1556                                   iCol = nextSlack[iCol];
1557                              }
1558                         }
1559                         //printf("now what?\n");
1560                         int n = 0;
1561                         bool basic = false;
1562                         if (lo > 0.0) {
1563                              // Add in positive
1564                              iCol = posSlack[i];
1565                              while (iCol >= 0) {
1566                                   double value = element[columnStart[iCol]];
1567                                   ratio[n] = cost[iCol] / value;
1568                                   which[n++] = iCol;
1569                                   iCol = nextSlack[iCol];
1570                              }
1571                              CoinSort_2(ratio, ratio + n, which);
1572                              for (int i = 0; i < n; i++) {
1573                                   iCol = which[i];
1574                                   double value = element[columnStart[iCol]];
1575                                   if (lo >= upper[iCol]*value) {
1576                                        value *= upper[iCol];
1577                                        sum += value;
1578                                        lo -= value;
1579                                        colsol[iCol] = upper[iCol];
1580                                   } else {
1581                                        value = lo / value;
1582                                        sum += lo;
1583                                        lo = 0.0;
1584                                        colsol[iCol] = value;
1585                                        model_->setColumnStatus(iCol, ClpSimplex::basic);
1586                                        basic = true;
1587                                   }
1588                                   if (lo < 1.0e-7)
1589                                        break;
1590                              }
1591                         } else if (up < 0.0) {
1592                              // Use lo so coding is more similar
1593                              lo = -up;
1594                              // Add in negative
1595                              iCol = negSlack[i];
1596                              while (iCol >= 0) {
1597                                   double value = -element[columnStart[iCol]];
1598                                   ratio[n] = cost[iCol] / value;
1599                                   which[n++] = iCol;
1600                                   iCol = nextSlack[iCol];
1601                              }
1602                              CoinSort_2(ratio, ratio + n, which);
1603                              for (int i = 0; i < n; i++) {
1604                                   iCol = which[i];
1605                                   double value = -element[columnStart[iCol]];
1606                                   if (lo >= upper[iCol]*value) {
1607                                        value *= upper[iCol];
1608                                        sum += value;
1609                                        lo -= value;
1610                                        colsol[iCol] = upper[iCol];
1611                                   } else {
1612                                        value = lo / value;
1613                                        sum += lo;
1614                                        lo = 0.0;
1615                                        colsol[iCol] = value;
1616                                        model_->setColumnStatus(iCol, ClpSimplex::basic);
1617                                        basic = true;
1618                                   }
1619                                   if (lo < 1.0e-7)
1620                                        break;
1621                              }
1622                         }
1623                         if (0) {
1624                              double sum2 = 0.0;
1625                              double x = 0.0;
1626                              for (int k = 0; k < ncols; k++) {
1627                                   CoinBigIndex j;
1628                                   double value = colsol[k];
1629                                   x += value * cost[k];
1630                                   for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
1631                                        int irow = row[j];
1632                                        if (irow == i)
1633                                             sum2 += element[j] * value;
1634                                   }
1635                              }
1636                              printf("after sum %g <= %g <= %g cost %.18g (sum = %g)\n",
1637                                     rowlower[i], sum2, rowupper[i], x, sum);
1638                         }
1639                         rowsol[i] = sum;
1640                         if (basic) {
1641                              if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1642                                   model_->setRowStatus(i, ClpSimplex::atLowerBound);
1643                              else
1644                                   model_->setRowStatus(i, ClpSimplex::atUpperBound);
1645                         }
1646                    } else {
1647                         int n = 0;
1648                         int iCol;
1649                         iCol = posSlack[i];
1650                         while (iCol >= 0) {
1651                              if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1652                                        colsol[iCol] < upper[iCol] - 1.0e-8) {
1653                                   model_->setColumnStatus(iCol, ClpSimplex::basic);
1654                                   n++;
1655                              }
1656                              iCol = nextSlack[iCol];
1657                         }
1658                         iCol = negSlack[i];
1659                         while (iCol >= 0) {
1660                              if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1661                                        colsol[iCol] < upper[iCol] - 1.0e-8) {
1662                                   model_->setColumnStatus(iCol, ClpSimplex::basic);
1663                                   n++;
1664                              }
1665                              iCol = nextSlack[iCol];
1666                         }
1667                         if (n) {
1668                              if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1669                                   model_->setRowStatus(i, ClpSimplex::atLowerBound);
1670                              else
1671                                   model_->setRowStatus(i, ClpSimplex::atUpperBound);
1672#ifdef CLP_INVESTIGATE
1673                              if (n > 1)
1674                                   printf("%d basic on row %d!\n", n, i);
1675#endif
1676                         }
1677                    }
1678               }
1679          }
1680          delete [] ratio;
1681          delete [] which;
1682#else
1683          for (i = 0; i < nrows; i++) {
1684               int n = 0;
1685               int iCol;
1686               iCol = posSlack[i];
1687               while (iCol >= 0) {
1688                    if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1689                              colsol[iCol] < upper[iCol] - 1.0e-8) {
1690                         model_->setColumnStatus(iCol, ClpSimplex::basic);
1691                         n++;
1692                    }
1693                    iCol = nextSlack[iCol];
1694               }
1695               iCol = negSlack[i];
1696               while (iCol >= 0) {
1697                    if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1698                              colsol[iCol] < upper[iCol] - 1.0e-8) {
1699                         model_->setColumnStatus(iCol, ClpSimplex::basic);
1700                         n++;
1701                    }
1702                    iCol = nextSlack[iCol];
1703               }
1704               if (n) {
1705                    if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1706                         model_->setRowStatus(i, ClpSimplex::atLowerBound);
1707                    else
1708                         model_->setRowStatus(i, ClpSimplex::atUpperBound);
1709#ifdef CLP_INVESTIGATE
1710                    if (n > 1)
1711                         printf("%d basic on row %d!\n", n, i);
1712#endif
1713               }
1714          }
1715#endif
1716     }
1717     double maxmin;
1718     if (model_->getObjSense() == -1.0) {
1719          maxmin = -1.0;
1720     } else {
1721          maxmin = 1.0;
1722     }
1723     bool justValuesPass = majorIterations_ > 1000000;
1724     if (slackStart >= 0) {
1725          for (i = 0; i < nrows; i++) {
1726               model_->setRowStatus(i, ClpSimplex::superBasic);
1727          }
1728          for (i = slackStart; i < slackEnd; i++) {
1729               model_->setColumnStatus(i, ClpSimplex::basic);
1730          }
1731     } else {
1732          /* still try and put singletons rather than artificials in basis */
1733          for (i = 0; i < nrows; i++) {
1734              model_->setRowStatus(i, ClpSimplex::basic);
1735          }
1736          int ninbas = 0;
1737          for (i = 0; i < ncols; i++) {
1738              if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) {
1739                    CoinBigIndex j = columnStart[i];
1740                    double value = element[j];
1741                    int irow = row[j];
1742                    double rlo = rowlower[irow];
1743                    double rup = rowupper[irow];
1744                    double clo = lower[i];
1745                    double cup = upper[i];
1746                    double csol = colsol[i];
1747                    /* adjust towards feasibility */
1748                    double move = 0.0;
1749                    if (rowsol[irow] > rup) {
1750                         move = (rup - rowsol[irow]) / value;
1751                         if (value > 0.0) {
1752                              /* reduce */
1753                              if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1754                         } else {
1755                              /* increase */
1756                              if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1757                         }
1758                    } else if (rowsol[irow] < rlo) {
1759                         move = (rlo - rowsol[irow]) / value;
1760                         if (value > 0.0) {
1761                              /* increase */
1762                              if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1763                         } else {
1764                              /* reduce */
1765                              if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1766                         }
1767                    } else {
1768                         /* move to improve objective */
1769                         if (cost[i]*maxmin > 0.0) {
1770                              if (value > 0.0) {
1771                                   move = (rlo - rowsol[irow]) / value;
1772                                   /* reduce */
1773                                   if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1774                              } else {
1775                                   move = (rup - rowsol[irow]) / value;
1776                                   /* increase */
1777                                   if (csol + move > cup) 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) move = CoinMax(0.0, cup - csol);
1784                              } else {
1785                                   move = (rlo - rowsol[irow]) / value;
1786                                   /* reduce */
1787                                   if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1788                              }
1789                         }
1790                    }
1791                    rowsol[irow] += move * value;
1792                    colsol[i] += move;
1793                    /* put in basis if row was artificial */
1794                    if (rup - rlo < 1.0e-7 && model_->getRowStatus(irow) == ClpSimplex::basic) {
1795                         model_->setRowStatus(irow, ClpSimplex::superBasic);
1796                         model_->setColumnStatus(i, ClpSimplex::basic);
1797                         ninbas++;
1798                    }
1799               }
1800            }
1801            /*printf("%d in basis\n",ninbas);*/
1802     }
1803     bool wantVector = false;
1804     if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
1805          // See if original wanted vector
1806          ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
1807          wantVector = clpMatrixO->wantsSpecialColumnCopy();
1808     }
1809     if ((strategy_&32768)!=0)
1810       allowInfeasible=true;
1811     if ((strategy_&65536)!=0)
1812       justValuesPass=true;
1813     //double * saveBounds=NULL;
1814     if (addAll < 3) {
1815          ClpPresolve pinfo;
1816          if (presolve) {
1817               double * rhs = new double[nrows];
1818               double * saveBounds=new double[2*ncols];
1819               char line[200];
1820               memcpy(saveBounds,lower,ncols*sizeof(double));
1821               memcpy(saveBounds+ncols,upper,ncols*sizeof(double));
1822               if (allowInfeasible) {
1823                    // fix up so will be feasible
1824                    const double * dual = model_->dualRowSolution();
1825                    for (i = 0; i < nrows; i++) 
1826                      rhs[i]=fabs(dual[i]);
1827                    std::sort(rhs,rhs+nrows);
1828                    int nSmall=nrows;
1829                    int nMedium=nrows;
1830                    double largest=rhs[nrows-1];
1831                    double small=CoinMax(1.0e-4,1.0e-5*largest);
1832                    small = CoinMin(small,1.0e-2);
1833                    double medium=small*100.0;
1834                    double * rowupper = model_->rowUpper();
1835                    double * rowlower = model_->rowLower();
1836                    // if tiny then drop row??
1837                    for (i = 0; i < nrows; i++) {
1838                      if (rhs[i]>=small) {
1839                        nSmall=i-1;
1840                        break;
1841                      }
1842                    }
1843                    for (; i < nrows; i++) {
1844                      if (rhs[i]>=medium) {
1845                        nMedium=i-1;
1846                        break;
1847                      }
1848                    }
1849                    printf("%d < %g, %d < %g, %d <= %g\n",
1850                           nSmall,small,nMedium-nSmall,medium,nrows-nMedium,largest);
1851                    memset(rhs, 0, nrows * sizeof(double));
1852                    int nFixed=0;
1853                    for (i = 0; i < ncols; i++) {
1854                      if (colsol[i]<lower[i]+1.0e-8) {
1855                        upper[i]=lower[i];
1856                        colsol[i]=lower[i];
1857                        nFixed++;
1858                      } else if (colsol[i]>upper[i]-1.0e-8) {
1859                        lower[i]=upper[i];
1860                        colsol[i]=lower[i];
1861                        nFixed++;
1862                      }
1863                    }
1864                    model_->clpMatrix()->times(1.0, colsol, rhs);
1865                    saveRowUpper = CoinCopyOfArray(rowupper, nrows);
1866                    saveRowLower = CoinCopyOfArray(rowlower, nrows);
1867                    double sum = 0.0;
1868                    for (i = 0; i < nrows; i++) {
1869                         if (rhs[i] > rowupper[i]) {
1870                              sum += rhs[i] - rowupper[i];
1871                         }
1872                         if (rhs[i] < rowlower[i]) {
1873                              sum += rowlower[i] - rhs[i];
1874                         }
1875                    }
1876                    double averageInfeasibility=sum/nrows;
1877                    double check=CoinMin(1.0e-3,0.1*averageInfeasibility);
1878                    int nFixedRows=0;
1879                    int nFreed=0;
1880#define MESS_UP 0
1881                    for (i = 0; i < nrows; i++) {
1882                         if (rowupper[i]>rowlower[i]+check) {
1883                           // look at distance and sign of dual
1884                           if (dual[i]<-medium&&rowupper[i]-rhs[i]<check) {
1885                             rowupper[i]=rhs[i];
1886                             rowlower[i]=rowupper[i];
1887                             nFixedRows++;
1888                           } else if (dual[i]>medium&&rhs[i]-rowlower[i]<check) {
1889                             rowlower[i]=rhs[i];
1890                             rowupper[i]=rowlower[i];
1891                             nFixedRows++;
1892                           } else if (fabs(dual[i])<small&&
1893                                      rhs[i]-rowlower[i]>check&&
1894                                      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                    double averageInfeasibility=sum/nrows;
1944                    sprintf(line,"sum of infeasibilities %g - average %g, %d fixed columns", 
1945                           sum,averageInfeasibility,nFixed);
1946               }
1947               const CoinMessages * messages = model_->messagesPointer();
1948               model_->messageHandler()->message(CLP_GENERAL, *messages)
1949                 << line
1950                 << CoinMessageEol;
1951               delete [] rhs;
1952               saveModel = model_;
1953               pinfo.setPresolveActions(pinfo.presolveActions() | 16384);
1954               model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
1955               if (saveBounds) {
1956                 memcpy(saveModel->columnLower(),saveBounds,ncols*sizeof(double));
1957                 memcpy(saveModel->columnUpper(),saveBounds+ncols,ncols*sizeof(double));
1958                 delete [] saveBounds;
1959               }
1960               if (model_&&(strategy_&262144)!=0) {
1961                 int nrows = model_->getNumRows();
1962                 int ncols = model_->getNumCols();
1963                 double * lower = model_->columnLower();
1964                 double * upper = model_->columnUpper();
1965                 const double * rowlower = model_->getRowLower();
1966                 const double * rowupper = model_->getRowUpper();
1967                 double * rowsol = model_->primalRowSolution();
1968                 double * colsol = model_->primalColumnSolution();;
1969                 int ninbas = 0;
1970                 int * which = new int[2*ncols+nrows];
1971                 double * dj = model_->dualColumnSolution();
1972                 for (int i=0;i<ncols;i++) {
1973                   dj[i]=-CoinMin(upper[i]-colsol[i],colsol[i]-lower[i]);
1974                   which[i]=i;
1975                 }
1976                 CoinSort_2(dj,dj+ncols,which);
1977                 ninbas=CoinMin(ncols,nrows);
1978                 int * columnIsBasic=which+ncols;
1979                 int * rowIsBasic=columnIsBasic+ncols;
1980                 for (int i=0;i<nrows+ncols;i++)
1981                   columnIsBasic[i]=-1;
1982                 for (int i=0;i<ninbas;i++) {
1983                   int iColumn=which[i];
1984                   columnIsBasic[iColumn]=i;
1985                 }
1986                 // factorize
1987                 CoinFactorization factor;
1988                 factor.pivotTolerance(0.1);
1989                 factor.setDenseThreshold(0);
1990                 int status=-1;
1991                 // If initial is too dense - then all-slack may be better
1992                 double areaFactor=1.0; // was 4.0
1993                 const CoinPackedMatrix * matrix = model_->matrix();
1994                 while (status) {
1995                   status=factor.factorize(*matrix,rowIsBasic,columnIsBasic,areaFactor);
1996                   if (status==-99) {
1997                     // put all slacks in
1998                     for (int i=0;i<nrows;i++) 
1999                       rowIsBasic[i]=i;
2000                     for (int i=0;i<ncols;i++) 
2001                      columnIsBasic[i]=-1;
2002                     break;
2003                   } else if (status==-1) {
2004                     factor.pivotTolerance(0.99);
2005                     // put all slacks in
2006                     for (int i=0;i<nrows;i++) 
2007                       rowIsBasic[i]=i;
2008                     for (int i=0;i<ncols;i++) {
2009                       int iRow=columnIsBasic[i];
2010                       if (iRow>=0)
2011                         rowIsBasic[iRow]=-1; // out
2012                     }
2013                   }
2014                 }
2015                 delete [] which;
2016                 for (int i = 0; i < nrows; i++) {
2017                   if (rowIsBasic[i]>=0) {
2018                     model_->setRowStatus(i, ClpSimplex::basic);
2019                   } else if (rowlower[i]==rowupper[i]) {
2020                     model_->setRowStatus(i, ClpSimplex::isFixed);
2021                   } else if (rowsol[i]-rowlower[i]<rowupper[i]-rowsol[i]) {
2022                     model_->setRowStatus(i, ClpSimplex::atLowerBound);
2023                   } else {
2024                     model_->setRowStatus(i, ClpSimplex::atUpperBound);
2025                   }
2026                 }
2027                 for (int i=0;i<ncols;i++) {
2028                   if (colsol[i]>upper[i]-1.0e-7||
2029                       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               }
2038          }
2039          if (model_) {
2040               // See if we want to go all way
2041               int oldSize=2*saveModel->numberRows()+saveModel->numberColumns();
2042               int newSize=2*model_->numberRows()+model_->numberColumns();
2043               if (oldSize*2>newSize*3)
2044                 justValuesPass=false;
2045               if (!wantVector) {
2046                    //#define TWO_GOES
2047#ifdef ABC_INHERIT
2048#ifndef TWO_GOES
2049                    model_->dealWithAbc(1,justValuesPass ? 3 : 1);
2050#else
2051                    model_->dealWithAbc(1,1 + 11);
2052#endif
2053#else
2054#ifndef TWO_GOES
2055                    model_->primal(justValuesPass ? 2 : 1);
2056#else
2057                    model_->primal(1 + 11);
2058#endif
2059#endif
2060               } else {
2061                    ClpMatrixBase * matrix = model_->clpMatrix();
2062                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
2063                    assert (clpMatrix);
2064                    clpMatrix->makeSpecialColumnCopy();
2065#ifdef ABC_INHERIT
2066                    model_->dealWithAbc(1,1);
2067#else
2068                    model_->primal(1);
2069#endif
2070                    clpMatrix->releaseSpecialColumnCopy();
2071               }
2072               if (presolve) {
2073                 if (!justValuesPass)
2074                   model_->primal(1);
2075                    pinfo.postsolve(true);
2076                    delete model_;
2077                    model_ = saveModel;
2078                    saveModel = NULL;
2079               } 
2080          } else {
2081               // not feasible
2082               addAll = 1;
2083               presolve = 0;
2084               model_ = saveModel;
2085               saveModel = NULL;
2086               if (justValuesPass)
2087#ifdef ABC_INHERIT
2088                    model_->dealWithAbc(1,3);
2089#else
2090                    model_->primal(2);
2091#endif
2092          }
2093          if (allowInfeasible) {
2094               CoinMemcpyN(saveRowUpper, nrows, model_->rowUpper());
2095               CoinMemcpyN(saveRowLower, nrows, model_->rowLower());
2096               delete [] saveRowUpper;
2097               delete [] saveRowLower;
2098               saveRowUpper = NULL;
2099               saveRowLower = NULL;
2100          }
2101          if (addAll < 2) {
2102               n = 0;
2103               if (!addAll ) {
2104                    /* could do scans to get a good number */
2105                    iteration = 1;
2106                    for (i = ordStart; i < ordEnd; i++) {
2107                         if (whenUsed[i] >= iteration) {
2108                              if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
2109                                   n++;
2110                                   upper[i] = saveUpper[i];
2111                                   lower[i] = saveLower[i];
2112                              }
2113                         }
2114                    }
2115               } else {
2116                    for (i = ordStart; i < ordEnd; i++) {
2117                         if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
2118                              n++;
2119                              upper[i] = saveUpper[i];
2120                              lower[i] = saveLower[i];
2121                         }
2122                    }
2123                    delete [] saveUpper;
2124                    delete [] saveLower;
2125                    saveUpper = NULL;
2126                    saveLower = NULL;
2127               }
2128#ifdef COIN_DEVELOP
2129               printf("Time so far %g, %d now added from previous iterations\n",
2130                      CoinCpuTime() - startTime, n);
2131#endif
2132               if (justValuesPass)
2133                    return;
2134               if (addAll)
2135                    presolve = 0;
2136               if (presolve) {
2137                    saveModel = model_;
2138                    model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
2139               } else {
2140                    presolve = 0;
2141               }
2142               if (!wantVector) {
2143#ifdef ABC_INHERIT
2144                    model_->dealWithAbc(1,1);
2145#else
2146                    model_->primal(1);
2147#endif
2148               } else {
2149                    ClpMatrixBase * matrix = model_->clpMatrix();
2150                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
2151                    assert (clpMatrix);
2152                    clpMatrix->makeSpecialColumnCopy();
2153#ifdef ABC_INHERIT
2154                    model_->dealWithAbc(1,1);
2155#else
2156                    model_->primal(1);
2157#endif
2158                    clpMatrix->releaseSpecialColumnCopy();
2159               }
2160               if (presolve) {
2161                    pinfo.postsolve(true);
2162                    delete model_;
2163                    model_ = saveModel;
2164                    saveModel = NULL;
2165               }
2166               if (!addAll) {
2167                    n = 0;
2168                    for (i = ordStart; i < ordEnd; i++) {
2169                         if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
2170                              n++;
2171                              upper[i] = saveUpper[i];
2172                              lower[i] = saveLower[i];
2173                         }
2174                    }
2175                    delete [] saveUpper;
2176                    delete [] saveLower;
2177                    saveUpper = NULL;
2178                    saveLower = NULL;
2179#ifdef COIN_DEVELOP
2180                    printf("Time so far %g, %d now added from previous iterations\n",
2181                           CoinCpuTime() - startTime, n);
2182#endif
2183               }
2184               if (presolve) {
2185                    saveModel = model_;
2186                    model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
2187               } else {
2188                    presolve = 0;
2189               }
2190               if (!wantVector) {
2191#ifdef ABC_INHERIT
2192                    model_->dealWithAbc(1,1);
2193#else
2194                    model_->primal(1);
2195#endif
2196               } else {
2197                    ClpMatrixBase * matrix = model_->clpMatrix();
2198                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
2199                    assert (clpMatrix);
2200                    clpMatrix->makeSpecialColumnCopy();
2201#ifdef ABC_INHERIT
2202                    model_->dealWithAbc(1,1);
2203#else
2204                    model_->primal(1);
2205#endif
2206                    clpMatrix->releaseSpecialColumnCopy();
2207               }
2208               if (presolve) {
2209                    pinfo.postsolve(true);
2210                    delete model_;
2211                    model_ = saveModel;
2212                    saveModel = NULL;
2213               }
2214          }
2215#ifdef COIN_DEVELOP
2216          printf("Total time in crossover %g\n", CoinCpuTime() - startTime);
2217#endif
2218          delete [] saveUpper;
2219          delete [] saveLower;
2220     }
2221#ifdef FEB_TRY
2222     model_->setSpecialOptions(saveOptions);
2223     model_->setPerturbation(savePerturbation);
2224#endif
2225     return ;
2226}
2227#endif
2228/*****************************************************************************/
2229
2230// Default contructor
2231Idiot::Idiot()
2232{
2233     model_ = NULL;
2234     maxBigIts_ = 3;
2235     maxIts_ = 5;
2236     logLevel_ = 1;
2237     logFreq_ = 100;
2238     maxIts2_ = 100;
2239     djTolerance_ = 1e-1;
2240     mu_ = 1e-4;
2241     drop_ = 5.0;
2242     exitDrop_ = -1.0e20;
2243     muFactor_ = 0.3333;
2244     stopMu_ = 1e-12;
2245     smallInfeas_ = 1e-1;
2246     reasonableInfeas_ = 1e2;
2247     muAtExit_ = 1.0e31;
2248     strategy_ = 8;
2249     lambdaIterations_ = 0;
2250     checkFrequency_ = 100;
2251     whenUsed_ = NULL;
2252     majorIterations_ = 30;
2253     exitFeasibility_ = -1.0;
2254     dropEnoughFeasibility_ = 0.02;
2255     dropEnoughWeighted_ = 0.01;
2256     // adjust
2257     double nrows = 10000.0;
2258     int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
2259     baseIts = baseIts / 10;
2260     baseIts *= 10;
2261     maxIts2_ = 200 + baseIts + 5;
2262     maxIts2_ = 100;
2263     reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
2264     lightWeight_ = 0;
2265}
2266// Constructor from model
2267Idiot::Idiot(OsiSolverInterface &model)
2268{
2269     model_ = & model;
2270     maxBigIts_ = 3;
2271     maxIts_ = 5;
2272     logLevel_ = 1;
2273     logFreq_ = 100;
2274     maxIts2_ = 100;
2275     djTolerance_ = 1e-1;
2276     mu_ = 1e-4;
2277     drop_ = 5.0;
2278     exitDrop_ = -1.0e20;
2279     muFactor_ = 0.3333;
2280     stopMu_ = 1e-12;
2281     smallInfeas_ = 1e-1;
2282     reasonableInfeas_ = 1e2;
2283     muAtExit_ = 1.0e31;
2284     strategy_ = 8;
2285     lambdaIterations_ = 0;
2286     checkFrequency_ = 100;
2287     whenUsed_ = NULL;
2288     majorIterations_ = 30;
2289     exitFeasibility_ = -1.0;
2290     dropEnoughFeasibility_ = 0.02;
2291     dropEnoughWeighted_ = 0.01;
2292     // adjust
2293     double nrows;
2294     if (model_)
2295          nrows = model_->getNumRows();
2296     else
2297          nrows = 10000.0;
2298     int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
2299     baseIts = baseIts / 10;
2300     baseIts *= 10;
2301     maxIts2_ = 200 + baseIts + 5;
2302     maxIts2_ = 100;
2303     reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
2304     lightWeight_ = 0;
2305}
2306// Copy constructor.
2307Idiot::Idiot(const Idiot &rhs)
2308{
2309     model_ = rhs.model_;
2310     if (model_ && rhs.whenUsed_) {
2311          int numberColumns = model_->getNumCols();
2312          whenUsed_ = new int [numberColumns];
2313          CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
2314     } else {
2315          whenUsed_ = NULL;
2316     }
2317     djTolerance_ = rhs.djTolerance_;
2318     mu_ = rhs.mu_;
2319     drop_ = rhs.drop_;
2320     muFactor_ = rhs.muFactor_;
2321     stopMu_ = rhs.stopMu_;
2322     smallInfeas_ = rhs.smallInfeas_;
2323     reasonableInfeas_ = rhs.reasonableInfeas_;
2324     exitDrop_ = rhs.exitDrop_;
2325     muAtExit_ = rhs.muAtExit_;
2326     exitFeasibility_ = rhs.exitFeasibility_;
2327     dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
2328     dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
2329     maxBigIts_ = rhs.maxBigIts_;
2330     maxIts_ = rhs.maxIts_;
2331     majorIterations_ = rhs.majorIterations_;
2332     logLevel_ = rhs.logLevel_;
2333     logFreq_ = rhs.logFreq_;
2334     checkFrequency_ = rhs.checkFrequency_;
2335     lambdaIterations_ = rhs.lambdaIterations_;
2336     maxIts2_ = rhs.maxIts2_;
2337     strategy_ = rhs.strategy_;
2338     lightWeight_ = rhs.lightWeight_;
2339}
2340// Assignment operator. This copies the data
2341Idiot &
2342Idiot::operator=(const Idiot & rhs)
2343{
2344     if (this != &rhs) {
2345          delete [] whenUsed_;
2346          model_ = rhs.model_;
2347          if (model_ && rhs.whenUsed_) {
2348               int numberColumns = model_->getNumCols();
2349               whenUsed_ = new int [numberColumns];
2350               CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
2351          } else {
2352               whenUsed_ = NULL;
2353          }
2354          djTolerance_ = rhs.djTolerance_;
2355          mu_ = rhs.mu_;
2356          drop_ = rhs.drop_;
2357          muFactor_ = rhs.muFactor_;
2358          stopMu_ = rhs.stopMu_;
2359          smallInfeas_ = rhs.smallInfeas_;
2360          reasonableInfeas_ = rhs.reasonableInfeas_;
2361          exitDrop_ = rhs.exitDrop_;
2362          muAtExit_ = rhs.muAtExit_;
2363          exitFeasibility_ = rhs.exitFeasibility_;
2364          dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
2365          dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
2366          maxBigIts_ = rhs.maxBigIts_;
2367          maxIts_ = rhs.maxIts_;
2368          majorIterations_ = rhs.majorIterations_;
2369          logLevel_ = rhs.logLevel_;
2370          logFreq_ = rhs.logFreq_;
2371          checkFrequency_ = rhs.checkFrequency_;
2372          lambdaIterations_ = rhs.lambdaIterations_;
2373          maxIts2_ = rhs.maxIts2_;
2374          strategy_ = rhs.strategy_;
2375          lightWeight_ = rhs.lightWeight_;
2376     }
2377     return *this;
2378}
2379Idiot::~Idiot()
2380{
2381     delete [] whenUsed_;
2382}
Note: See TracBrowser for help on using the repository browser.