source: stable/BSP/Clp/src/Idiot.cpp @ 1251

Last change on this file since 1251 was 1251, checked in by stefan, 12 years ago

disable some printf's per default

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