source: trunk/Clp/src/Idiot.cpp @ 1115

Last change on this file since 1115 was 1115, checked in by forrest, 12 years ago

mainly fix isFixed bug when reusing factorization

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