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

Last change on this file since 1054 was 1054, checked in by forrest, 13 years ago

changes to compile under gcc 4.3

  • 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    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    printf("%d slacks\n",numberSlacks);
478    oddSlacks=true;
479    int extra = (int) (nrows*sizeof(double)/sizeof(int));
480    whenUsed_=new int[2*ncols+2*nrows+extra];
481    int * posSlack = whenUsed_+ncols;
482    int * negSlack = posSlack+nrows;
483    int * nextSlack = negSlack + nrows;
484    for (i=0;i<nrows;i++) {
485      posSlack[i]=-1;
486      negSlack[i]=-1;
487    }
488    for (i=0;i<ncols;i++) 
489      nextSlack[i]=-1;
490    for (i=0;i<ncols;i++) {
491      if (columnLength[i]==1) {
492        CoinBigIndex j=columnStart[i];
493        int iRow = row[j];
494        if (element[j]>0.0) {
495          if (posSlack[iRow]==-1) {
496            posSlack[iRow]=i;
497          } else {
498            int iCol = posSlack[iRow];
499            while (nextSlack[iCol]>=0)
500              iCol = nextSlack[iCol];
501            nextSlack[iCol]=i;
502          }
503        } else {
504          if (negSlack[iRow]==-1) {
505            negSlack[iRow]=i;
506          } else {
507            int iCol = negSlack[iRow];
508            while (nextSlack[iCol]>=0)
509              iCol = nextSlack[iCol];
510            nextSlack[iCol]=i;
511          }
512        }
513      }
514    }
515    // now sort
516    for (i=0;i<nrows;i++) {
517      int iCol;
518      iCol = posSlack[i];
519      if (iCol>=0) {
520        CoinBigIndex j=columnStart[iCol];
521#ifndef NDEBUG
522        int iRow = row[j];
523#endif
524        assert (element[j]>0.0);
525        assert (iRow==i);
526        dj[0]= cost[iCol]/element[j];
527        whenUsed_[0]=iCol;
528        int n=1;
529        while (nextSlack[iCol]>=0) {
530          iCol = nextSlack[iCol];
531          CoinBigIndex j=columnStart[iCol];
532#ifndef NDEBUG
533          int iRow = row[j];
534#endif
535          assert (element[j]>0.0);
536          assert (iRow==i);
537          dj[n]= cost[iCol]/element[j];
538          whenUsed_[n++]=iCol;
539        }
540        for (j=0;j<n;j++) {
541          int jCol = whenUsed_[j];
542          nextSlack[jCol]=-2;
543        }
544        CoinSort_2(dj,dj+n,whenUsed_);
545        // put back
546        iCol = whenUsed_[0];
547        posSlack[i]=iCol;
548        for (j=1;j<n;j++) {
549          int jCol = whenUsed_[j];
550          nextSlack[iCol]=jCol;
551          iCol=jCol;
552        }
553      }
554      iCol = negSlack[i];
555      if (iCol>=0) {
556        CoinBigIndex j=columnStart[iCol];
557#ifndef NDEBUG
558        int iRow = row[j];
559#endif
560        assert (element[j]<0.0);
561        assert (iRow==i);
562        dj[0]= -cost[iCol]/element[j];
563        whenUsed_[0]=iCol;
564        int n=1;
565        while (nextSlack[iCol]>=0) {
566          iCol = nextSlack[iCol];
567          CoinBigIndex j=columnStart[iCol];
568#ifndef NDEBUG
569          int iRow = row[j];
570#endif
571          assert (element[j]<0.0);
572          assert (iRow==i);
573          dj[n]= -cost[iCol]/element[j];
574          whenUsed_[n++]=iCol;
575        }
576        for (j=0;j<n;j++) {
577          int jCol = whenUsed_[j];
578          nextSlack[jCol]=-2;
579        }
580        CoinSort_2(dj,dj+n,whenUsed_);
581        // put back
582        iCol = whenUsed_[0];
583        negSlack[i]=iCol;
584        for (j=1;j<n;j++) {
585          int jCol = whenUsed_[j];
586          nextSlack[iCol]=jCol;
587          iCol=jCol;
588        }
589      }
590    }
591  }
592  whenUsed=whenUsed_;
593  if (model_->getObjSense()==-1.0) {
594    maxmin=-1.0;
595  } else {
596    maxmin=1.0;
597  }
598  model_->getDblParam(OsiObjOffset,offset);
599  if (!maxIts2) maxIts2=maxIts;
600  strategy=strategy_;
601  strategy &= 3;
602  memset(lambda,0,nrows*sizeof(double));
603  slackStart=countCostedSlacks(model_);
604  if (slackStart>=0) {
605    printf("This model has costed slacks\n");
606    slackEnd=slackStart+nrows;
607    if (slackStart) {
608      ordStart=0;
609      ordEnd=slackStart;
610    } else {
611      ordStart=nrows;
612      ordEnd=ncols;
613    }
614  } else {
615    slackEnd=slackStart;
616    ordStart=0;
617    ordEnd=ncols;
618  }
619  if (offset&&logLevel>2) {
620    printf("** Objective offset is %g\n",offset);
621  }
622  /* compute reasonable solution cost */
623  for (i=0;i<nrows;i++) {
624    rowsol[i]=1.0e31;
625  }
626  for (i=0;i<ncols;i++) {
627    CoinBigIndex j;
628    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
629      if (element[j]!=1.0) {
630        allOnes=0;
631        break;
632      }
633    }
634  }
635  if (allOnes) {
636    elemXX=NULL;
637  } else {
638    elemXX=element;
639  }
640  // Do scaling if wanted
641  bool scaled=false;
642#ifndef OSI_IDIOT
643  if ((strategy_&32)!=0&&!allOnes) {
644    if (model_->scalingFlag()>0)
645      scaled = model_->clpMatrix()->scale(model_)==0;
646    if (scaled) {
647      const double * rowScale = model_->rowScale();
648      const double * columnScale = model_->columnScale();
649      double * oldLower = lower;
650      double * oldUpper = upper;
651      double * oldCost = cost;
652      lower = new double[ncols];
653      upper = new double[ncols];
654      cost = new double[ncols];
655      memcpy(lower,oldLower,ncols*sizeof(double));
656      memcpy(upper,oldUpper,ncols*sizeof(double));
657      memcpy(cost,oldCost,ncols*sizeof(double));
658      int icol,irow;
659      for (icol=0;icol<ncols;icol++) {
660        double multiplier = 1.0/columnScale[icol];
661        if (lower[icol]>-1.0e50)
662          lower[icol] *= multiplier;
663        if (upper[icol]<1.0e50)
664          upper[icol] *= multiplier;
665        colsol[icol] *= multiplier;
666        cost[icol] *= columnScale[icol];
667      }
668      memcpy(rowlower,model_->rowLower(),nrows*sizeof(double));
669      for (irow=0;irow<nrows;irow++) {
670        double multiplier = rowScale[irow];
671        if (rowlower[irow]>-1.0e50)
672          rowlower[irow] *= multiplier;
673        if (rowupper[irow]<1.0e50)
674          rowupper[irow] *= multiplier;
675        rowsol[irow] *= multiplier;
676      }
677      int length = columnStart[ncols-1]+columnLength[ncols-1];
678      double * elemYY = new double[length];
679      for (i=0;i<ncols;i++) {
680        CoinBigIndex j;
681        double scale = columnScale[i];
682        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
683          int irow=row[j];
684          elemYY[j] = element[j]*scale*rowScale[irow];
685        }
686      }
687      elemXX=elemYY;
688    }
689  }
690#endif
691  for (i=0;i<ncols;i++) {
692    CoinBigIndex j;
693    double dd=columnLength[i];
694    dd=cost[i]/dd;
695    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
696      int irow=row[j];
697      if (dd<rowsol[irow]) {
698        rowsol[irow]=dd;
699      }
700    }
701  }
702  d2=0.0;
703  for (i=0;i<nrows;i++) {
704    d2+=rowsol[i];
705  }
706  d2*=2.0; /* for luck */
707 
708  d2=d2/((double) (4*nrows+8000));
709  d2*=0.5; /* halve with more flexible method */
710  if (d2<5.0) d2=5.0;
711  if (djExit==0.0) {
712    djExit=d2;
713  }
714  if ((saveStrategy&4)!=0) {
715    /* go to relative tolerances - first small */
716    djExit=1.0e-10;
717    djFlag=1.0e-5;
718    drop=1.0e-10;
719  }
720  memset(whenUsed,0,ncols*sizeof(int));
721  strategy=strategies[strategy];
722  if ((saveStrategy&8)!=0) strategy |= 64; /* don't allow large theta */
723  memcpy(saveSol,colsol,ncols*sizeof(double));
724 
725  lastResult=IdiSolve(nrows,ncols,rowsol ,colsol,pi,
726                       dj,cost,rowlower,rowupper,
727                       lower,upper,elemXX,row,columnStart,columnLength,lambda,
728                       0,mu,drop,
729                       maxmin,offset,strategy,djTol,djExit,djFlag);
730  // update whenUsed_
731  n = cleanIteration(iteration, ordStart,ordEnd,
732                     colsol,  lower,  upper,
733                     rowlower, rowupper,
734                     cost, elemXX, fixTolerance,lastResult.objval,lastResult.infeas);
735  if ((strategy_&16384)!=0) {
736    int * posSlack = whenUsed_+ncols;
737    int * negSlack = posSlack+nrows;
738    int * nextSlack = negSlack + nrows;
739    double * rowsol2 = (double *) (nextSlack+ncols);
740    for (i=0;i<nrows;i++)
741      rowsol[i] += rowsol2[i];
742  }
743  if ((logLevel_&1)!=0) {
744#ifndef OSI_IDIOT
745    if (!handler) {
746#endif
747      printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
748             iteration,lastResult.infeas,lastResult.objval,mu,lastResult.iteration,n);
749#ifndef OSI_IDIOT
750    } else {
751      handler->message(CLP_IDIOT_ITERATION,*messages)
752        <<iteration<<lastResult.infeas<<lastResult.objval<<mu<<lastResult.iteration<<n
753        <<CoinMessageEol;
754    }
755#endif
756  }
757  int numberBaseTrys=0; // for first time
758  int numberAway=-1;
759  iterationTotal = lastResult.iteration;
760  firstInfeas=lastResult.infeas;
761  if ((strategy_&1024)!=0) reasonableInfeas=0.5*firstInfeas;
762  if (lastResult.infeas<reasonableInfeas) lastResult.infeas=reasonableInfeas;
763  double keepinfeas=1.0e31;
764  double lastInfeas=1.0e31;
765  double bestInfeas=1.0e31;
766  while ((mu>stopMu&&lastResult.infeas>smallInfeas)||
767         (lastResult.infeas<=smallInfeas&&
768         dropping(lastResult,exitDrop,smallInfeas,&badIts))||
769         checkIteration<2||lambdaIteration<lambdaIterations_) {
770    if (lastResult.infeas<=exitFeasibility_)
771      break; 
772    iteration++;
773    checkIteration++;
774    if (lastResult.infeas<=smallInfeas&&lastResult.objval<bestFeasible) {
775      bestFeasible=lastResult.objval;
776    }
777    if (lastResult.infeas+mu*lastResult.objval<bestWeighted) {
778      bestWeighted=lastResult.objval+mu*lastResult.objval;
779    }
780    if ((saveStrategy&4096)) strategy |=256;
781    if ((saveStrategy&4)!=0&&iteration>2) {
782      /* go to relative tolerances */
783      double weighted=10.0+fabs(lastWeighted);
784      djExit=djTolerance_*weighted;
785      djFlag=2.0*djExit;
786      drop=drop_*weighted;
787      djTol=0.01*djExit;
788    }
789    memcpy(saveSol,colsol,ncols*sizeof(double));
790    result=IdiSolve(nrows,ncols,rowsol ,colsol,pi,dj,
791                     cost,rowlower,rowupper,
792                     lower,upper,elemXX,row,columnStart,columnLength,lambda,
793                     maxIts,mu,drop,
794                     maxmin,offset,strategy,djTol,djExit,djFlag);
795    n = cleanIteration(iteration, ordStart,ordEnd,
796                       colsol,  lower,  upper,
797                       rowlower, rowupper,
798                       cost, elemXX, fixTolerance,result.objval,result.infeas);
799    if ((strategy_&16384)!=0) {
800      int * posSlack = whenUsed_+ncols;
801      int * negSlack = posSlack+nrows;
802      int * nextSlack = negSlack + nrows;
803      double * rowsol2 = (double *) (nextSlack+ncols);
804      for (i=0;i<nrows;i++)
805        rowsol[i] += rowsol2[i];
806    }
807    if ((logLevel_&1)!=0) {
808#ifndef OSI_IDIOT
809      if (!handler) {
810#endif
811        printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
812               iteration,result.infeas,result.objval,mu,result.iteration,n);
813#ifndef OSI_IDIOT
814      } else {
815        handler->message(CLP_IDIOT_ITERATION,*messages)
816          <<iteration<<result.infeas<<result.objval<<mu<<result.iteration<<n
817          <<CoinMessageEol;
818      }
819#endif
820    }
821    if (iteration>50&&n==numberAway&&result.infeas<1.0e-4) {
822      printf("infeas small %g\n",result.infeas);
823      break; // not much happening
824    }
825    if (lightWeight_==1&&iteration>10&&result.infeas>1.0&&maxIts!=7) {
826      if (lastInfeas!=bestInfeas&&CoinMin(result.infeas,lastInfeas)>0.95*bestInfeas)
827        majorIterations_ = CoinMin(majorIterations_,iteration); // not getting feasible
828    }
829    lastInfeas = result.infeas;
830    numberAway=n;
831    keepinfeas = result.infeas;
832    lastWeighted=result.weighted;
833    iterationTotal += result.iteration;
834    if (iteration==1) {
835      if ((strategy_&1024)!=0&&mu<1.0e-10) 
836        result.infeas=firstInfeas*0.8;
837      if (majorIterations_>=50||dropEnoughFeasibility_<=0.0)
838        result.infeas *= 0.8;
839      if (result.infeas>firstInfeas*0.9
840          &&result.infeas>reasonableInfeas) {
841        iteration--;
842        if (majorIterations_<50)
843          mu*=1.0e-1;
844        else
845          mu*=0.7;
846        bestFeasible=1.0e31;
847        bestWeighted=1.0e60;
848        numberBaseTrys++;
849        if (mu<1.0e-30||(numberBaseTrys>10&&lightWeight_)) {
850          // back to all slack basis
851          lightWeight_=2;
852          break;
853        }
854        memcpy(colsol,saveSol,ncols*sizeof(double));
855      } else {
856        maxIts=maxIts2;
857        checkIteration=0;
858        if ((strategy_&1024)!=0) mu *= 1.0e-1;
859      }
860    } else {
861    }
862    bestInfeas=CoinMin(bestInfeas,result.infeas);
863    if (iteration) {
864      /* this code is in to force it to terminate sometime */
865      double changeMu=factor;
866      if ((saveStrategy&64)!=0) {
867        keepinfeas=0.0; /* switch off ranga's increase */
868        fakeSmall=smallInfeas;
869      } else {
870        fakeSmall=-1.0;
871      }
872      saveLambdaScale=0.0;
873      if (result.infeas>reasonableInfeas||
874          (nTry+1==maxBigIts&&result.infeas>fakeSmall)) {
875        if (result.infeas>lastResult.infeas*(1.0-dropEnoughFeasibility_)||
876            nTry+1==maxBigIts||
877            (result.infeas>lastResult.infeas*0.9
878             &&result.weighted>lastResult.weighted
879             -dropEnoughWeighted_*CoinMax(fabs(lastResult.weighted),fabs(result.weighted)))) {
880          mu*=changeMu;
881          if ((saveStrategy&32)!=0&&result.infeas<reasonableInfeas&&0) {
882            reasonableInfeas=CoinMax(smallInfeas,reasonableInfeas*sqrt(changeMu));
883            printf("reasonable infeas now %g\n",reasonableInfeas);
884          }
885          result.weighted=1.0e60;
886          nTry=0;
887          bestFeasible=1.0e31;
888          bestWeighted=1.0e60;
889          checkIteration=0;
890          lambdaIteration=0;
891#define LAMBDA
892#ifdef LAMBDA
893          if ((saveStrategy&2048)==0) {
894            memset(lambda,0,nrows*sizeof(double));
895          }
896#else
897          memset(lambda,0,nrows*sizeof(double));
898#endif
899        } else {
900          nTry++;
901        }
902      } else if (lambdaIterations_>=0) {
903        /* update lambda  */
904        double scale=1.0/mu;
905        int i,nnz=0;
906        saveLambdaScale=scale;
907         lambdaIteration++;
908         if ((saveStrategy&4)==0) drop = drop_/50.0;
909         if (lambdaIteration>4 && 
910            (((lambdaIteration%10)==0 && smallInfeas<keepinfeas) ||
911             ((lambdaIteration%5)==0 && 1.5*smallInfeas<keepinfeas))) {
912           //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
913           smallInfeas *= 1.5;
914         }
915         if ((saveStrategy&2048)==0) {
916           for (i=0;i<nrows;i++) {
917             if (lambda[i]) nnz++;
918             lambda[i]+= scale*rowsol[i];
919           }
920         } else {
921           nnz=1;
922#ifdef LAMBDA
923           for (i=0;i<nrows;i++) {
924             lambda[i]+= scale*rowsol[i];
925           }
926#else
927           for (i=0;i<nrows;i++) {
928             lambda[i] = scale*rowsol[i];
929           }
930           for (i=0;i<ncols;i++) {
931             CoinBigIndex j;
932             double value=cost[i]*maxmin;
933             for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
934               int irow=row[j];
935               value+=element[j]*lambda[irow];
936             }
937             cost[i]=value*maxmin;
938           }
939           for (i=0;i<nrows;i++) {
940             offset+=lambda[i]*rowupper[i];
941             lambda[i]=0.0;
942           }
943#ifdef DEBUG
944           printf("offset %g\n",offset);
945#endif
946           model_->setDblParam(OsiObjOffset,offset);
947#endif
948         }
949        nTry++;
950        if (!nnz) {
951          bestFeasible=1.0e32;
952          bestWeighted=1.0e60;
953          checkIteration=0;
954          result.weighted=1.0e31;
955        }
956#ifdef DEBUG
957        double trueCost=0.0;
958        for (i=0;i<ncols;i++) {
959          int j;
960          trueCost+=cost[i]*colsol[i];
961        }
962        printf("True objective %g\n",trueCost-offset);
963#endif
964      } else {
965        nTry++;
966      }
967      lastResult=result;
968      if (result.infeas<reasonableInfeas&&!belowReasonable) {
969        belowReasonable=1;
970        bestFeasible=1.0e32;
971        bestWeighted=1.0e60;
972        checkIteration=0;
973        result.weighted=1.0e31;
974      }
975    }
976    if (iteration>=majorIterations_) {
977      // If not feasible and crash then dive dive dive
978      if (mu>1.0e-12&&result.infeas>1.0&&majorIterations_<40) {
979        mu=1.0e-30;
980        majorIterations_=iteration+1;
981        stopMu=0.0;
982      } else {
983        if (logLevel>2) 
984          printf("Exiting due to number of major iterations\n");
985        break;
986      }
987    }
988  }
989  majorIterations_ = saveMajorIterations;
990#ifndef OSI_IDIOT
991  if (scaled) {
992    // Scale solution and free arrays
993    const double * rowScale = model_->rowScale();
994    const double * columnScale = model_->columnScale();
995    int icol,irow;
996    for (icol=0;icol<ncols;icol++) {
997      colsol[icol] *= columnScale[icol];
998      saveSol[icol] *= columnScale[icol];
999      dj[icol] /= columnScale[icol];
1000    }
1001    for (irow=0;irow<nrows;irow++) {
1002      rowsol[irow] /= rowScale[irow];
1003      pi[irow] *= rowScale[irow];
1004    }
1005    // Don't know why getting Microsoft problems
1006#if defined (_MSC_VER)
1007    delete [] ( double *) elemXX;
1008#else
1009    delete [] elemXX;
1010#endif
1011    model_->setRowScale(NULL);
1012    model_->setColumnScale(NULL);
1013    delete [] lower;
1014    delete [] upper;
1015    delete [] cost;
1016    lower = model_->columnLower();
1017    upper = model_->columnUpper();
1018    cost = model_->objective();
1019    //rowlower = model_->rowLower();
1020  }
1021#endif
1022#define TRYTHIS
1023#ifdef TRYTHIS
1024  if ((saveStrategy&2048)!=0) {
1025    double offset;
1026    model_->getDblParam(OsiObjOffset,offset);
1027    for (i=0;i<ncols;i++) {
1028      CoinBigIndex j;
1029      double djval=cost[i]*maxmin;
1030      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1031        int irow=row[j];
1032        djval -= element[j]*lambda[irow];
1033      }
1034      cost[i]=djval;
1035    }
1036    for (i=0;i<nrows;i++) {
1037      offset+=lambda[i]*rowupper[i];
1038    }
1039    model_->setDblParam(OsiObjOffset,offset);
1040  }
1041#endif
1042  if (saveLambdaScale) {
1043    /* back off last update */
1044    for (i=0;i<nrows;i++) {
1045      lambda[i]-= saveLambdaScale*rowsol[i];
1046    }
1047  }
1048  muAtExit_=mu;
1049  // For last iteration make as feasible as possible
1050  if (oddSlacks)
1051    strategy_ |= 16384;
1052  // not scaled
1053  n = cleanIteration(iteration, ordStart,ordEnd,
1054                     colsol,  lower,  upper,
1055                     model_->rowLower(), model_->rowUpper(),
1056                     cost, element, fixTolerance,lastResult.objval,lastResult.infeas);
1057#if 0
1058  if ((logLevel&1)==0||(strategy_&16384)!=0) {
1059    printf(
1060            "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
1061            iteration,mu,lastResult.infeas,lastResult.objval,n);
1062  }
1063#endif
1064#ifndef OSI_IDIOT
1065  model_->setSumPrimalInfeasibilities(lastResult.infeas);
1066#endif
1067  // Put back more feasible solution
1068  double saveInfeas[]={0.0,0.0};
1069  for (int iSol=0;iSol<3;iSol++) {
1070    const double * solution = iSol ? colsol : saveSol;
1071    if (iSol==2&&saveInfeas[0]<saveInfeas[1]) {
1072      // put back best solution
1073      memcpy(colsol,saveSol,ncols*sizeof(double));
1074    }
1075    double large=0.0;
1076    int i;
1077    memset(rowsol,0,nrows*sizeof(double));
1078    for (i=0;i<ncols;i++) {
1079      CoinBigIndex j;
1080      double value=solution[i];
1081      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1082        int irow=row[j];
1083        rowsol[irow] += element[j]*value;
1084      }
1085    }
1086    for (i=0;i<nrows;i++) {
1087      if (rowsol[i] > rowupper[i]) {
1088        double diff=rowsol[i]-rowupper[i];
1089        if (diff>large) 
1090          large=diff;
1091      } else if (rowsol[i] < rowlower[i]) {
1092        double diff=rowlower[i]-rowsol[i];
1093        if (diff>large) 
1094          large=diff;
1095      } 
1096    }
1097    if (iSol<2)
1098      saveInfeas[iSol]=large;
1099    if (logLevel>2)
1100      printf("largest infeasibility is %g\n",large);
1101  }
1102  /* subtract out lambda */
1103  for (i=0;i<nrows;i++) {
1104    pi[i]-=lambda[i];
1105  }
1106  for (i=0;i<ncols;i++) {
1107    CoinBigIndex j;
1108    double djval=cost[i]*maxmin;
1109    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1110      int irow=row[j];
1111      djval -= element[j]*pi[irow];
1112    }
1113    dj[i]=djval;
1114  }
1115  if ((strategy_&1024)!=0) {
1116    double ratio = ((double) ncols)/((double) nrows);
1117    printf("col/row ratio %g infeas ratio %g\n",ratio,lastResult.infeas/firstInfeas);
1118    if (lastResult.infeas>0.01*firstInfeas*ratio) {
1119      strategy_ &= (~1024);
1120      printf(" - layer off\n");
1121    } else {
1122      printf(" - layer on\n");
1123    }
1124  }
1125  delete [] saveSol;
1126  delete [] lambda;
1127  // save solution
1128  // duals not much use - but save anyway
1129#ifndef OSI_IDIOT
1130  memcpy(model_->primalRowSolution(),rowsol,nrows*sizeof(double));
1131  memcpy(model_->primalColumnSolution(),colsol,ncols*sizeof(double));
1132  memcpy(model_->dualRowSolution(),pi,nrows*sizeof(double));
1133  memcpy(model_->dualColumnSolution(),dj,ncols*sizeof(double));
1134#else
1135  model_->setColSolution(colsol);
1136  model_->setRowPrice(pi);
1137  delete [] cost;
1138#endif
1139  delete [] rowsol;
1140  delete [] colsol;
1141  delete [] pi;
1142  delete [] dj;
1143  delete [] rowlower;
1144  delete [] rowupper;
1145  return ;
1146}
1147#ifndef OSI_IDIOT
1148void
1149Idiot::crossOver(int mode)
1150{
1151  if (lightWeight_==2) {
1152    // total failure
1153    model_->allSlackBasis();
1154    return;
1155  }
1156  double fixTolerance=IDIOT_FIX_TOLERANCE;
1157  double startTime = CoinCpuTime();
1158  ClpSimplex * saveModel=NULL;
1159  ClpMatrixBase * matrix = model_->clpMatrix();
1160  const int * row = matrix->getIndices();
1161  const CoinBigIndex * columnStart = matrix->getVectorStarts();
1162  const int * columnLength = matrix->getVectorLengths(); 
1163  const double * element = matrix->getElements();
1164  const double * rowupper = model_->getRowUpper();
1165  int nrows=model_->getNumRows();
1166  int ncols=model_->getNumCols();
1167  double * rowsol, * colsol;
1168  // different for Osi
1169  double * lower = model_->columnLower();
1170  double * upper = model_->columnUpper();
1171  const double * rowlower= model_->getRowLower();
1172  int * whenUsed=whenUsed_;
1173  rowsol= model_->primalRowSolution();
1174  colsol= model_->primalColumnSolution();;
1175  double * cost=model_->objective();
1176
1177  int slackEnd,ordStart,ordEnd;
1178  int slackStart = countCostedSlacks(model_);
1179
1180  int addAll = mode&7;
1181  int presolve=0;
1182
1183  double djTolerance = djTolerance_;
1184  if (djTolerance>0.0&&djTolerance<1.0)
1185    djTolerance=1.0;
1186  int iteration;
1187  int i, n=0;
1188  double ratio=1.0;
1189  double objValue=0.0;
1190  if ((strategy_&128)!=0) {
1191    fixTolerance=SMALL_IDIOT_FIX_TOLERANCE;
1192  }
1193  if ((mode&16)!=0&&addAll<3) presolve=1;
1194  double * saveUpper = NULL;
1195  double * saveLower = NULL;
1196  double * saveRowUpper = NULL;
1197  double * saveRowLower = NULL;
1198  bool allowInfeasible = (strategy_&8192)!=0;
1199  if (addAll<3) {
1200    saveUpper = new double [ncols];
1201    saveLower = new double [ncols];
1202    memcpy(saveUpper,upper,ncols*sizeof(double));
1203    memcpy(saveLower,lower,ncols*sizeof(double));
1204    if (allowInfeasible) {
1205      saveRowUpper = new double [nrows];
1206      saveRowLower = new double [nrows];
1207      memcpy(saveRowUpper,rowupper,nrows*sizeof(double));
1208      memcpy(saveRowLower,rowlower,nrows*sizeof(double));
1209      double averageInfeas = model_->sumPrimalInfeasibilities()/((double) model_->numberRows());
1210      fixTolerance = CoinMax(fixTolerance,1.0e-5*averageInfeas);
1211    }
1212  }
1213  if (slackStart>=0) {
1214    slackEnd=slackStart+nrows;
1215    if (slackStart) {
1216      ordStart=0;
1217      ordEnd=slackStart;
1218    } else {
1219      ordStart=nrows;
1220      ordEnd=ncols;
1221    }
1222  } else {
1223    slackEnd=slackStart;
1224    ordStart=0;
1225    ordEnd=ncols;
1226  }
1227  /* get correct rowsol (without known slacks) */
1228  memset(rowsol,0,nrows*sizeof(double));
1229  for (i=ordStart;i<ordEnd;i++) {
1230    CoinBigIndex j;
1231    double value=colsol[i];
1232    if (value<lower[i]+fixTolerance) {
1233      value=lower[i];
1234      colsol[i]=value;
1235    }
1236    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1237      int irow=row[j];
1238      rowsol[irow]+=value*element[j];
1239    }
1240  }
1241  if (slackStart>=0) {
1242    for (i=0;i<nrows;i++) {
1243      if (ratio*rowsol[i]>rowlower[i]&&rowsol[i]>1.0e-8) {
1244        ratio=rowlower[i]/rowsol[i];
1245      }
1246    }
1247    for (i=0;i<nrows;i++) {
1248      rowsol[i]*=ratio;
1249    }
1250    for (i=ordStart;i<ordEnd;i++) {
1251      double value=colsol[i]*ratio;
1252      colsol[i]=value;
1253      objValue+=value*cost[i];
1254    }
1255    for (i=0;i<nrows;i++) {
1256      double value=rowlower[i]-rowsol[i];
1257      colsol[i+slackStart]=value;
1258      objValue+=value*cost[i+slackStart];
1259    }
1260    printf("New objective after scaling %g\n",objValue);
1261  }
1262#if 0
1263   maybe put back - but just get feasible ?
1264  // If not many fixed then just exit
1265  int numberFixed=0;
1266  for (i=ordStart;i<ordEnd;i++) {
1267    if (colsol[i]<lower[i]+fixTolerance)
1268      numberFixed++;
1269    else if (colsol[i]>upper[i]-fixTolerance)
1270      numberFixed++;
1271  }
1272  if (numberFixed<ncols/2) {
1273    addAll=3;
1274    presolve=0;
1275  }
1276#endif
1277  model_->createStatus();
1278  /* addAll
1279     0 - chosen,all used, all
1280     1 - chosen, all
1281     2 - all
1282     3 - do not do anything  - maybe basis
1283  */
1284  for (i=ordStart;i<ordEnd;i++) {
1285    if (addAll<2) {
1286      if (colsol[i]<lower[i]+fixTolerance) {
1287        upper[i]=lower[i];
1288        colsol[i]=lower[i];
1289      } else if (colsol[i]>upper[i]-fixTolerance) {
1290        lower[i]=upper[i];
1291        colsol[i]=upper[i];
1292      }
1293    }
1294    model_->setColumnStatus(i,ClpSimplex::superBasic);
1295  }
1296  if ((strategy_&16384)!=0) {
1297    // put in basis
1298    int * posSlack = whenUsed_+ncols;
1299    int * negSlack = posSlack+nrows;
1300    int * nextSlack = negSlack + nrows;
1301    for (i=0;i<nrows;i++) {
1302      int n=0;
1303      int iCol;
1304      iCol =posSlack[i];
1305      if (iCol>=0) {
1306        if (colsol[iCol]>lower[iCol]+1.0e-8&&
1307            colsol[iCol]<upper[iCol]+1.0e-8) {
1308          model_->setColumnStatus(iCol,ClpSimplex::basic);
1309          n++;
1310        }
1311        while (nextSlack[iCol]>=0) {
1312          iCol = nextSlack[iCol];
1313          if (colsol[iCol]>lower[iCol]+1.0e-8&&
1314              colsol[iCol]<upper[iCol]+1.0e-8) {
1315            model_->setColumnStatus(iCol,ClpSimplex::basic);
1316            n++;
1317          }
1318        }
1319      }
1320      iCol =negSlack[i];
1321      if (iCol>=0) {
1322        if (colsol[iCol]>lower[iCol]+1.0e-8&&
1323            colsol[iCol]<upper[iCol]+1.0e-8) {
1324          model_->setColumnStatus(iCol,ClpSimplex::basic);
1325          n++;
1326        }
1327        while (nextSlack[iCol]>=0) {
1328          iCol = nextSlack[iCol];
1329          if (colsol[iCol]>lower[iCol]+1.0e-8&&
1330              colsol[iCol]<upper[iCol]+1.0e-8) {
1331            model_->setColumnStatus(iCol,ClpSimplex::basic);
1332            n++;
1333          }
1334        }
1335      }
1336      if (n) {
1337        if (fabs(rowsol[i]-rowlower[i])<fabs(rowsol[i]-rowupper[i]))
1338          model_->setRowStatus(i,ClpSimplex::atLowerBound);
1339        else
1340          model_->setRowStatus(i,ClpSimplex::atUpperBound);
1341        if (n>1)
1342          printf("%d basic on row %d!\n",n,i);
1343      }
1344    }
1345  }
1346  double maxmin;
1347  if (model_->getObjSense()==-1.0) {
1348    maxmin=-1.0;
1349  } else {
1350    maxmin=1.0;
1351  }
1352  if (slackStart>=0) {
1353    for (i=0;i<nrows;i++) {
1354      model_->setRowStatus(i,ClpSimplex::superBasic);
1355    }
1356    for (i=slackStart;i<slackEnd;i++) {
1357      model_->setColumnStatus(i,ClpSimplex::basic);
1358    }
1359  } else {
1360    /* still try and put singletons rather than artificials in basis */
1361    int ninbas=0;
1362    for (i=0;i<nrows;i++) {
1363      model_->setRowStatus(i,ClpSimplex::basic);
1364    }
1365    for (i=0;i<ncols;i++) {
1366      if (columnLength[i]==1&&upper[i]>lower[i]+1.0e-5) {
1367        CoinBigIndex j =columnStart[i];
1368        double value=element[j];
1369        int irow=row[j];
1370        double rlo=rowlower[irow];
1371        double rup=rowupper[irow];
1372        double clo=lower[i];
1373        double cup=upper[i];
1374        double csol=colsol[i];
1375        /* adjust towards feasibility */
1376        double move=0.0;
1377        if (rowsol[irow]>rup) {
1378          move=(rup-rowsol[irow])/value;
1379          if (value>0.0) {
1380            /* reduce */
1381            if (csol+move<clo) move=CoinMin(0.0,clo-csol);
1382          } else {
1383            /* increase */
1384            if (csol+move>cup) move=CoinMax(0.0,cup-csol);
1385          }
1386        } else if (rowsol[irow]<rlo) {
1387          move=(rlo-rowsol[irow])/value;
1388          if (value>0.0) {
1389            /* increase */
1390            if (csol+move>cup) move=CoinMax(0.0,cup-csol);
1391          } else {
1392            /* reduce */
1393            if (csol+move<clo) move=CoinMin(0.0,clo-csol);
1394          }
1395        } else {
1396          /* move to improve objective */
1397          if (cost[i]*maxmin>0.0) {
1398            if (value>0.0) {
1399              move=(rlo-rowsol[irow])/value;
1400              /* reduce */
1401              if (csol+move<clo) move=CoinMin(0.0,clo-csol);
1402            } else {
1403              move=(rup-rowsol[irow])/value;
1404              /* increase */
1405              if (csol+move>cup) move=CoinMax(0.0,cup-csol);
1406            }
1407          } else if (cost[i]*maxmin<0.0) {
1408            if (value>0.0) {
1409              move=(rup-rowsol[irow])/value;
1410              /* increase */
1411              if (csol+move>cup) move=CoinMax(0.0,cup-csol);
1412            } else {
1413              move=(rlo-rowsol[irow])/value;
1414              /* reduce */
1415              if (csol+move<clo) move=CoinMin(0.0,clo-csol);
1416            }
1417          }
1418        }
1419        rowsol[irow] +=move*value;
1420        colsol[i]+=move;
1421        /* put in basis if row was artificial */
1422        if (rup-rlo<1.0e-7&&model_->getRowStatus(irow)==ClpSimplex::basic) {
1423          model_->setRowStatus(irow,ClpSimplex::superBasic);
1424          model_->setColumnStatus(i,ClpSimplex::basic);
1425          ninbas++;
1426        }
1427      }
1428    }
1429    /*printf("%d in basis\n",ninbas);*/
1430  }
1431  bool wantVector=false;
1432  if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
1433    // See if original wanted vector
1434    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
1435    wantVector = clpMatrixO->wantsSpecialColumnCopy();
1436  }
1437  if (addAll<3) {
1438    ClpPresolve pinfo;
1439    if (presolve) {
1440      if (allowInfeasible) {
1441        // fix up so will be feasible
1442        double * rhs = new double[nrows];
1443        memset(rhs,0,nrows*sizeof(double));
1444        model_->clpMatrix()->times(1.0,colsol,rhs);
1445        double * rowupper = model_->rowUpper();
1446        double * rowlower= model_->rowLower();
1447        double sum = 0.0;
1448        for (i=0;i<nrows;i++) {
1449          if (rhs[i]>rowupper[i]) {
1450            sum += rhs[i]-rowupper[i];
1451            rowupper[i]=rhs[i];
1452          }
1453          if (rhs[i]<rowlower[i]) {
1454            sum += rowlower[i]-rhs[i];
1455            rowlower[i]=rhs[i];
1456          }
1457        }
1458        printf("sum of infeasibilities %g\n",sum);
1459        delete [] rhs;
1460      }
1461      saveModel = model_;
1462      pinfo.setPresolveActions(pinfo.presolveActions()|16384);
1463      model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1464    }
1465    if (model_) {
1466      if (!wantVector) {
1467        model_->primal(1);
1468      } else {
1469        ClpMatrixBase * matrix = model_->clpMatrix();
1470        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1471        assert (clpMatrix);
1472        clpMatrix->makeSpecialColumnCopy();
1473        model_->primal(1);
1474        clpMatrix->releaseSpecialColumnCopy();
1475      }
1476      if (presolve) {
1477        pinfo.postsolve(true);
1478        delete model_;
1479        model_ = saveModel;
1480        saveModel=NULL;
1481      }
1482    } else {
1483      // not feasible
1484      addAll=1;
1485      presolve=0;
1486      model_ = saveModel;
1487      saveModel=NULL;
1488    }
1489    if (allowInfeasible) {
1490      memcpy(model_->rowUpper(),saveRowUpper,nrows*sizeof(double));
1491      memcpy(model_->rowLower(),saveRowLower,nrows*sizeof(double));
1492      delete [] saveRowUpper;
1493      delete [] saveRowLower;
1494      saveRowUpper = NULL;
1495      saveRowLower = NULL;
1496    }
1497    if (addAll<2) {
1498      n=0;
1499      if (!addAll ) {
1500        /* could do scans to get a good number */
1501        iteration=1;
1502        for (i=ordStart;i<ordEnd;i++) {
1503          if (whenUsed[i]>=iteration) {
1504            if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1505              n++;
1506              upper[i]=saveUpper[i];
1507              lower[i]=saveLower[i];
1508            }
1509          }
1510        }
1511      } else {
1512        for (i=ordStart;i<ordEnd;i++) {
1513          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1514            n++;
1515            upper[i]=saveUpper[i];
1516            lower[i]=saveLower[i];
1517          }
1518        }
1519        delete [] saveUpper;
1520        delete [] saveLower;
1521        saveUpper=NULL;
1522        saveLower=NULL;
1523      }
1524      printf("Time so far %g, %d now added from previous iterations\n",
1525             CoinCpuTime()-startTime,n);
1526      if (addAll)
1527        presolve=0;
1528      if (presolve) {
1529        saveModel = model_;
1530        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1531      } else {
1532        presolve=0;
1533      }
1534      if (!wantVector) {
1535        model_->primal(1);
1536      } else {
1537        ClpMatrixBase * matrix = model_->clpMatrix();
1538        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1539        assert (clpMatrix);
1540        clpMatrix->makeSpecialColumnCopy();
1541        model_->primal(1);
1542        clpMatrix->releaseSpecialColumnCopy();
1543      }
1544      if (presolve) {
1545        pinfo.postsolve(true);
1546        delete model_;
1547        model_ = saveModel;
1548        saveModel=NULL;
1549      }
1550      if (!addAll) {
1551        n=0;
1552        for (i=ordStart;i<ordEnd;i++) {
1553          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1554            n++;
1555            upper[i]=saveUpper[i];
1556            lower[i]=saveLower[i];
1557          }
1558        }
1559        delete [] saveUpper;
1560        delete [] saveLower;
1561        saveUpper=NULL;
1562        saveLower=NULL;
1563        printf("Time so far %g, %d now added from previous iterations\n",
1564               CoinCpuTime()-startTime,n);
1565      }
1566      if (presolve) {
1567        saveModel = model_;
1568        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1569      } else {
1570        presolve=0;
1571      }
1572      if (!wantVector) {
1573        model_->primal(1);
1574      } else {
1575        ClpMatrixBase * matrix = model_->clpMatrix();
1576        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1577        assert (clpMatrix);
1578        clpMatrix->makeSpecialColumnCopy();
1579        model_->primal(1);
1580        clpMatrix->releaseSpecialColumnCopy();
1581      }
1582      if (presolve) {
1583        pinfo.postsolve(true);
1584        delete model_;
1585        model_ = saveModel;
1586        saveModel=NULL;
1587      }
1588    }
1589    printf("Total time in crossover %g\n", CoinCpuTime()-startTime);
1590    delete [] saveUpper;
1591    delete [] saveLower;
1592  }
1593  return ;
1594}
1595#endif
1596/*****************************************************************************/
1597
1598// Default contructor
1599Idiot::Idiot()
1600{
1601  model_ = NULL;
1602  maxBigIts_ = 3;
1603  maxIts_ = 5;
1604  logLevel_ = 1; 
1605  logFreq_ = 100;
1606  maxIts2_ = 100;
1607  djTolerance_ = 1e-1;
1608  mu_ = 1e-4;
1609  drop_ = 5.0;
1610  exitDrop_=-1.0e20;
1611  muFactor_ = 0.3333;
1612  stopMu_ = 1e-12;
1613  smallInfeas_ = 1e-1;
1614  reasonableInfeas_ = 1e2;
1615  muAtExit_ =1.0e31;
1616  strategy_ =8;
1617  lambdaIterations_ =0;
1618  checkFrequency_ =100;
1619  whenUsed_ = NULL;
1620  majorIterations_ =30;
1621  exitFeasibility_ =-1.0;
1622  dropEnoughFeasibility_ =0.02;
1623  dropEnoughWeighted_ =0.01;
1624  // adjust
1625  double nrows=10000.0;
1626  int baseIts =(int) sqrt((double)nrows);
1627  baseIts =baseIts/10;
1628  baseIts *= 10;
1629  maxIts2_ =200+baseIts+5;
1630  maxIts2_=100;
1631  reasonableInfeas_ =((double) nrows)*0.05;
1632  lightWeight_=0;
1633}
1634// Constructor from model
1635Idiot::Idiot(OsiSolverInterface &model)
1636{
1637  model_ = & model;
1638  maxBigIts_ = 3;
1639  maxIts_ = 5;
1640  logLevel_ = 1; 
1641  logFreq_ = 100;
1642  maxIts2_ = 100;
1643  djTolerance_ = 1e-1;
1644  mu_ = 1e-4;
1645  drop_ = 5.0;
1646  exitDrop_=-1.0e20;
1647  muFactor_ = 0.3333;
1648  stopMu_ = 1e-12;
1649  smallInfeas_ = 1e-1;
1650  reasonableInfeas_ = 1e2;
1651  muAtExit_ =1.0e31;
1652  strategy_ =8;
1653  lambdaIterations_ =0;
1654  checkFrequency_ =100;
1655  whenUsed_ = NULL;
1656  majorIterations_ =30;
1657  exitFeasibility_ =-1.0;
1658  dropEnoughFeasibility_ =0.02;
1659  dropEnoughWeighted_ =0.01;
1660  // adjust
1661  double nrows;
1662  if (model_)
1663    nrows=model_->getNumRows();
1664  else
1665    nrows=10000.0;
1666  int baseIts =(int) sqrt((double)nrows);
1667  baseIts =baseIts/10;
1668  baseIts *= 10;
1669  maxIts2_ =200+baseIts+5;
1670  maxIts2_=100;
1671  reasonableInfeas_ =((double) nrows)*0.05;
1672  lightWeight_=0;
1673}
1674// Copy constructor.
1675Idiot::Idiot(const Idiot &rhs)
1676{
1677  model_ = rhs.model_;
1678  if (model_&&rhs.whenUsed_) {
1679    int numberColumns = model_->getNumCols();
1680    whenUsed_ = new int [numberColumns];
1681    memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1682  } else {
1683    whenUsed_=NULL;
1684  }
1685  djTolerance_ = rhs.djTolerance_;
1686  mu_ = rhs.mu_;
1687  drop_ = rhs.drop_;
1688  muFactor_ = rhs.muFactor_;
1689  stopMu_ = rhs.stopMu_;
1690  smallInfeas_ = rhs.smallInfeas_;
1691  reasonableInfeas_ = rhs.reasonableInfeas_;
1692  exitDrop_ = rhs.exitDrop_;
1693  muAtExit_ = rhs.muAtExit_;
1694  exitFeasibility_ = rhs.exitFeasibility_;
1695  dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1696  dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1697  maxBigIts_ = rhs.maxBigIts_;
1698  maxIts_ = rhs.maxIts_;
1699  majorIterations_ = rhs.majorIterations_;
1700  logLevel_ = rhs.logLevel_;
1701  logFreq_ = rhs.logFreq_;
1702  checkFrequency_ = rhs.checkFrequency_;
1703  lambdaIterations_ = rhs.lambdaIterations_;
1704  maxIts2_ = rhs.maxIts2_;
1705  strategy_ = rhs.strategy_;
1706  lightWeight_=rhs.lightWeight_;
1707}
1708// Assignment operator. This copies the data
1709Idiot & 
1710Idiot::operator=(const Idiot & rhs)
1711{
1712  if (this != &rhs) {
1713    delete [] whenUsed_;
1714    model_ = rhs.model_;
1715    if (model_&&rhs.whenUsed_) {
1716      int numberColumns = model_->getNumCols();
1717      whenUsed_ = new int [numberColumns];
1718      memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1719    } else {
1720      whenUsed_=NULL;
1721    }
1722    djTolerance_ = rhs.djTolerance_;
1723    mu_ = rhs.mu_;
1724    drop_ = rhs.drop_;
1725    muFactor_ = rhs.muFactor_;
1726    stopMu_ = rhs.stopMu_;
1727    smallInfeas_ = rhs.smallInfeas_;
1728    reasonableInfeas_ = rhs.reasonableInfeas_;
1729    exitDrop_ = rhs.exitDrop_;
1730    muAtExit_ = rhs.muAtExit_;
1731    exitFeasibility_ = rhs.exitFeasibility_;
1732    dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1733    dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1734    maxBigIts_ = rhs.maxBigIts_;
1735    maxIts_ = rhs.maxIts_;
1736    majorIterations_ = rhs.majorIterations_;
1737    logLevel_ = rhs.logLevel_;
1738    logFreq_ = rhs.logFreq_;
1739    checkFrequency_ = rhs.checkFrequency_;
1740    lambdaIterations_ = rhs.lambdaIterations_;
1741    maxIts2_ = rhs.maxIts2_;
1742    strategy_ = rhs.strategy_;
1743    lightWeight_=rhs.lightWeight_;
1744  }
1745  return *this;
1746}
1747Idiot::~Idiot()
1748{
1749  delete [] whenUsed_;
1750}
Note: See TracBrowser for help on using the repository browser.