source: branches/devel/Clp/src/Idiot.cpp @ 989

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

this may be a mistake

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