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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.1 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include <stdio.h>
6#include <stdarg.h>
7#include <stdlib.h>
8#include <math.h>
9#include "ClpPresolve.hpp"
10#include "Idiot.hpp"
11#include "CoinTime.hpp"
12#include "CoinMessageHandler.hpp"
13#include "CoinHelperFunctions.hpp"
14// Redefine stuff for Clp
15#ifndef OSI_IDIOT
16#include "ClpMessage.hpp"
17#define OsiObjOffset ClpObjOffset
18#endif
19/**** strategy 4 - drop, exitDrop and djTolerance all relative:
20For first two major iterations these are small.  Then:
21
22drop - exit a major iteration if drop over 5*checkFrequency < this is
23used as info->drop*(10.0+fabs(last weighted objective))
24
25exitDrop - exit idiot if feasible and drop < this is
26used as info->exitDrop*(10.0+fabs(last objective))
27
28djExit - exit a major iteration if largest dj (averaged over 5 checks)
29drops below this - used as info->djTolerance*(10.0+fabs(last weighted objective)
30
31djFlag - mostly skip variables with bad dj worse than this => 2*djExit
32
33djTol - only look at variables with dj better than this => 0.01*djExit
34****************/
35
36#define IDIOT_FIX_TOLERANCE 1e-6
37#define SMALL_IDIOT_FIX_TOLERANCE 1e-10
38int 
39Idiot::dropping(IdiotResult result,
40                    double tolerance,
41                    double small,
42                    int *nbad)
43{
44  if (result.infeas<=small) {
45    double value=CoinMax(fabs(result.objval),fabs(result.dropThis))+1.0;
46    if (result.dropThis>tolerance*value) {
47      *nbad=0;
48      return 1;
49    } else {
50      (*nbad)++;
51      if (*nbad>4) {
52        return 0;
53      } else {
54        return 1;
55      }
56    }
57  } else {
58    *nbad=0;
59    return 1;
60  }
61}
62
63/* returns -1 or start of costed slacks */
64 static int countCostedSlacks(OsiSolverInterface * model)
65{
66#ifdef OSI_IDIOT
67  const CoinPackedMatrix * matrix = model->getMatrixByCol();
68#else
69  ClpMatrixBase * matrix = model->clpMatrix();
70#endif
71  const int * row = matrix->getIndices();
72  const CoinBigIndex * columnStart = matrix->getVectorStarts();
73  const int * columnLength = matrix->getVectorLengths(); 
74  const double * element = matrix->getElements();
75  const double * rowupper = model->getRowUpper();
76  int nrows=model->getNumRows();
77  int ncols=model->getNumCols();
78  int slackStart = ncols-nrows;
79  int nSlacks=nrows;
80  int i;
81 
82  if (ncols<=nrows) return -1;
83  while (1) {
84    for (i=0;i<nrows;i++) {
85      int j=i+slackStart;
86      CoinBigIndex k=columnStart[j];
87      if (columnLength[j]==1) {
88        if (row[k]!=i||element[k]!=1.0) {
89          nSlacks=0;
90          break;
91        }
92      } else {
93        nSlacks=0;
94        break;
95      }
96      if (rowupper[i]<=0.0) {
97        nSlacks=0;
98        break;
99      }
100    }
101    if (nSlacks||!slackStart) break;
102    slackStart=0;
103  }
104  if (!nSlacks) slackStart=-1;
105  return slackStart;
106}
107void
108Idiot::crash(int numberPass, CoinMessageHandler * handler,const CoinMessages *messages)
109{
110  // lightweight options
111  int numberColumns = model_->getNumCols();
112  const double * objective = model_->getObjCoefficients();
113  int nnzero=0;
114  double sum=0.0;
115  int i;
116  for (i=0;i<numberColumns;i++) {
117    if (objective[i]) {
118      sum += fabs(objective[i]);
119      nnzero++;
120    }
121  }
122  sum /= (double) (nnzero+1);
123  maxIts_=2;
124  if (numberPass<=0)
125    // Cast to double to avoid VACPP complaining
126    majorIterations_=(int)(2+log10((double)(numberColumns+1)));
127  else
128    majorIterations_=numberPass;
129  // If mu not changed then compute
130  if (mu_==1e-4)
131    mu_= CoinMax(1.0e-3,sum*1.0e-5);
132  if (!lightWeight_) {
133    maxIts2_=105;
134  } else if (lightWeight_==1) {
135    mu_ *= 1000.0;
136    maxIts2_=23;
137  } else if (lightWeight_==2) {
138    maxIts2_=11;
139  } else {
140    maxIts2_=23;
141  }
142  //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
143  solve2(handler,messages);
144#ifndef OSI_IDIOT
145  double averageInfeas = model_->sumPrimalInfeasibilities()/((double) model_->numberRows());
146  if (averageInfeas<0.01&&(strategy_&512)!=0) 
147    crossOver(16+1); 
148  else
149    crossOver(3);
150#endif
151}
152void
153Idiot::solve()
154{
155  CoinMessages dummy;
156  solve2(NULL,&dummy);
157}
158void
159Idiot::solve2(CoinMessageHandler * handler,const CoinMessages * messages)
160{
161  int strategy=0;
162  double d2;
163  int i,n;
164  int allOnes=1;
165  int iteration=0;
166  int iterationTotal=0;
167  int nTry=0; /* number of tries at same weight */
168  double fixTolerance=IDIOT_FIX_TOLERANCE;
169  int maxBigIts=maxBigIts_;
170  int maxIts=maxIts_;
171  int logLevel=logLevel_;
172  int saveMajorIterations = majorIterations_;
173  if (handler) {
174    if (handler->logLevel()>0&&handler->logLevel()<3)
175      logLevel=1;
176    else if (!handler->logLevel())
177      logLevel=0;
178    else
179      logLevel=7;
180  }
181  double djExit=djTolerance_;
182  double djFlag=1.0+100.0*djExit;
183  double djTol=0.00001;
184  double mu =mu_;
185  double drop=drop_;
186  int maxIts2=maxIts2_;
187  double factor=muFactor_;
188  double smallInfeas=smallInfeas_;
189  double reasonableInfeas=reasonableInfeas_;
190  double stopMu=stopMu_;
191  double maxmin,offset;
192  double lastWeighted=1.0e50;
193  double exitDrop=exitDrop_;
194  double fakeSmall=smallInfeas;
195  double firstInfeas;
196  int badIts=0;
197  int slackStart,slackEnd,ordStart,ordEnd;
198  int checkIteration=0;
199  int lambdaIteration=0;
200  int belowReasonable=0; /* set if ever gone below reasonable infeas */
201  double bestWeighted=1.0e60;
202  double bestFeasible=1.0e60; /* best solution while feasible */
203  IdiotResult result,lastResult;
204  int saveStrategy=strategy_;
205  const int strategies[]={0,2,128};
206  double saveLambdaScale=0.0;
207  if ((saveStrategy&128)!=0) {
208    fixTolerance=SMALL_IDIOT_FIX_TOLERANCE;
209  }
210#ifdef OSI_IDIOT
211  const CoinPackedMatrix * matrix = model_->getMatrixByCol();
212#else
213  ClpMatrixBase * matrix = model_->clpMatrix();
214#endif
215  const int * row = matrix->getIndices();
216  const CoinBigIndex * columnStart = matrix->getVectorStarts();
217  const int * columnLength = matrix->getVectorLengths(); 
218  const double * element = matrix->getElements();
219  int nrows=model_->getNumRows();
220  int ncols=model_->getNumCols();
221  double * rowsol, * colsol;
222  double * pi, * dj;
223#ifndef OSI_IDIOT
224  double * cost = model_->objective();
225  double * lower = model_->columnLower();
226  double * upper = model_->columnUpper();
227#else
228  double * cost = new double [ncols];
229  memcpy( cost, model_->getObjCoefficients(), ncols*sizeof(double));
230  const double * lower = model_->getColLower();
231  const double * upper = model_->getColUpper();
232#endif
233  const double *elemXX;
234  double * saveSol;
235  double * rowupper= new double[nrows]; // not const as modified
236  memcpy(rowupper,model_->getRowUpper(),nrows*sizeof(double));
237  double * rowlower= new double[nrows]; // not const as modified
238  memcpy(rowlower,model_->getRowLower(),nrows*sizeof(double));
239  int * whenUsed;
240  double * lambda;
241  saveSol=new double[ncols];
242  lambda=new double [nrows];
243  rowsol= new double[nrows];
244  colsol= new double [ncols];
245  memcpy(colsol,model_->getColSolution(),ncols*sizeof(double));
246  pi= new double[nrows];
247  dj=new double[ncols];
248  delete [] whenUsed_;
249  whenUsed=whenUsed_=new int[ncols];
250  if (model_->getObjSense()==-1.0) {
251    maxmin=-1.0;
252  } else {
253    maxmin=1.0;
254  }
255  model_->getDblParam(OsiObjOffset,offset);
256  if (!maxIts2) maxIts2=maxIts;
257  strategy=strategy_;
258  strategy &= 3;
259  memset(lambda,0,nrows*sizeof(double));
260  slackStart=countCostedSlacks(model_);
261  if (slackStart>=0) {
262    printf("This model has costed slacks\n");
263    slackEnd=slackStart+nrows;
264    if (slackStart) {
265      ordStart=0;
266      ordEnd=slackStart;
267    } else {
268      ordStart=nrows;
269      ordEnd=ncols;
270    }
271  } else {
272    slackEnd=slackStart;
273    ordStart=0;
274    ordEnd=ncols;
275  }
276  if (offset&&logLevel>2) {
277    printf("** Objective offset is %g\n",offset);
278  }
279  /* compute reasonable solution cost */
280  for (i=0;i<nrows;i++) {
281    rowsol[i]=1.0e31;
282  }
283  for (i=0;i<ncols;i++) {
284    CoinBigIndex j;
285    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
286      if (element[j]!=1.0) {
287        allOnes=0;
288        break;
289      }
290    }
291  }
292  if (allOnes) {
293    elemXX=NULL;
294  } else {
295    elemXX=element;
296  }
297  // Do scaling if wanted
298  bool scaled=false;
299#ifndef OSI_IDIOT
300  if ((strategy_&32)!=0&&!allOnes) {
301    if (model_->scalingFlag()>0)
302      scaled = model_->clpMatrix()->scale(model_)==0;
303    if (scaled) {
304      const double * rowScale = model_->rowScale();
305      const double * columnScale = model_->columnScale();
306      double * oldLower = lower;
307      double * oldUpper = upper;
308      double * oldCost = cost;
309      lower = new double[ncols];
310      upper = new double[ncols];
311      cost = new double[ncols];
312      memcpy(lower,oldLower,ncols*sizeof(double));
313      memcpy(upper,oldUpper,ncols*sizeof(double));
314      memcpy(cost,oldCost,ncols*sizeof(double));
315      int icol,irow;
316      for (icol=0;icol<ncols;icol++) {
317        double multiplier = 1.0/columnScale[icol];
318        if (lower[icol]>-1.0e50)
319          lower[icol] *= multiplier;
320        if (upper[icol]<1.0e50)
321          upper[icol] *= multiplier;
322        colsol[icol] *= multiplier;
323        cost[icol] *= columnScale[icol];
324      }
325      memcpy(rowlower,model_->rowLower(),nrows*sizeof(double));
326      for (irow=0;irow<nrows;irow++) {
327        double multiplier = rowScale[irow];
328        if (rowlower[irow]>-1.0e50)
329          rowlower[irow] *= multiplier;
330        if (rowupper[irow]<1.0e50)
331          rowupper[irow] *= multiplier;
332        rowsol[irow] *= multiplier;
333      }
334      int length = columnStart[ncols-1]+columnLength[ncols-1];
335      double * elemYY = new double[length];
336      for (i=0;i<ncols;i++) {
337        CoinBigIndex j;
338        double scale = columnScale[i];
339        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
340          int irow=row[j];
341          elemYY[j] = element[j]*scale*rowScale[irow];
342        }
343      }
344      elemXX=elemYY;
345    }
346  }
347#endif
348  for (i=0;i<ncols;i++) {
349    CoinBigIndex j;
350    double dd=columnLength[i];
351    dd=cost[i]/dd;
352    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
353      int irow=row[j];
354      if (dd<rowsol[irow]) {
355        rowsol[irow]=dd;
356      }
357    }
358  }
359  d2=0.0;
360  for (i=0;i<nrows;i++) {
361    d2+=rowsol[i];
362  }
363  d2*=2.0; /* for luck */
364 
365  d2=d2/((double) (4*nrows+8000));
366  d2*=0.5; /* halve with more flexible method */
367  if (d2<5.0) d2=5.0;
368  if (djExit==0.0) {
369    djExit=d2;
370  }
371  if ((saveStrategy&4)!=0) {
372    /* go to relative tolerances - first small */
373    djExit=1.0e-10;
374    djFlag=1.0e-5;
375    drop=1.0e-10;
376  }
377  memset(whenUsed,0,ncols*sizeof(int));
378  strategy=strategies[strategy];
379  if ((saveStrategy&8)!=0) strategy |= 64; /* don't allow large theta */
380  memcpy(saveSol,colsol,ncols*sizeof(double));
381 
382  lastResult=IdiSolve(nrows,ncols,rowsol ,colsol,pi,
383                       dj,cost,rowlower,rowupper,
384                       lower,upper,elemXX,row,columnStart,columnLength,lambda,
385                       0,mu,drop,
386                       maxmin,offset,strategy,djTol,djExit,djFlag);
387  n=0;
388  for (i=ordStart;i<ordEnd;i++) {
389    if (colsol[i]>lower[i]+fixTolerance) {
390      if (colsol[i]<upper[i]-fixTolerance) {
391        n++;
392      } else {
393        colsol[i]=upper[i];
394      }
395      whenUsed[i]=iteration;
396    } else {
397      colsol[i]=lower[i];
398    }
399  }
400  if ((logLevel_&1)!=0) {
401#ifndef OSI_IDIOT
402    if (!handler) {
403#endif
404      printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
405             iteration,lastResult.infeas,lastResult.objval,mu,lastResult.iteration,n);
406#ifndef OSI_IDIOT
407    } else {
408      handler->message(CLP_IDIOT_ITERATION,*messages)
409        <<iteration<<lastResult.infeas<<lastResult.objval<<mu<<lastResult.iteration<<n
410        <<CoinMessageEol;
411    }
412#endif
413  }
414  int numberBaseTrys=0; // for first time
415  int numberAway=-1;
416  iterationTotal = lastResult.iteration;
417  firstInfeas=lastResult.infeas;
418  if ((strategy_&1024)!=0) reasonableInfeas=0.5*firstInfeas;
419  if (lastResult.infeas<reasonableInfeas) lastResult.infeas=reasonableInfeas;
420  double keepinfeas=1.0e31;
421  double lastInfeas=1.0e31;
422  double bestInfeas=1.0e31;
423  while ((mu>stopMu&&lastResult.infeas>smallInfeas)||
424         (lastResult.infeas<=smallInfeas&&
425         dropping(lastResult,exitDrop,smallInfeas,&badIts))||
426         checkIteration<2||lambdaIteration<lambdaIterations_) {
427    if (lastResult.infeas<=exitFeasibility_)
428      break; 
429    iteration++;
430    checkIteration++;
431    if (lastResult.infeas<=smallInfeas&&lastResult.objval<bestFeasible) {
432      bestFeasible=lastResult.objval;
433    }
434    if ((saveStrategy&4096)) strategy |=256;
435    if ((saveStrategy&4)!=0&&iteration>2) {
436      /* go to relative tolerances */
437      double weighted=10.0+fabs(lastWeighted);
438      djExit=djTolerance_*weighted;
439      djFlag=2.0*djExit;
440      drop=drop_*weighted;
441      djTol=0.01*djExit;
442    }
443    result=IdiSolve(nrows,ncols,rowsol ,colsol,pi,dj,
444                     cost,rowlower,rowupper,
445                     lower,upper,elemXX,row,columnStart,columnLength,lambda,
446                     maxIts,mu,drop,
447                     maxmin,offset,strategy,djTol,djExit,djFlag);
448    n=0;
449    for (i=ordStart;i<ordEnd;i++) {
450      if (colsol[i]>lower[i]+fixTolerance) {
451        if (colsol[i]<upper[i]-fixTolerance) {
452          n++;
453        } else {
454          colsol[i]=upper[i];
455        }
456        whenUsed[i]=iteration;
457      } else {
458        colsol[i]=lower[i];
459      }
460    }
461    if ((logLevel_&1)!=0) {
462#ifndef OSI_IDIOT
463      if (!handler) {
464#endif
465        printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
466               iteration,result.infeas,result.objval,mu,result.iteration,n);
467#ifndef OSI_IDIOT
468      } else {
469        handler->message(CLP_IDIOT_ITERATION,*messages)
470          <<iteration<<result.infeas<<result.objval<<mu<<result.iteration<<n
471          <<CoinMessageEol;
472      }
473#endif
474    }
475    if (iteration>50&&n==numberAway&&result.infeas<1.0e-4)
476      break; // not much happening
477    if (lightWeight_==1&&iteration>10&&result.infeas>1.0&&maxIts!=7) {
478      if (lastInfeas!=bestInfeas&&CoinMin(result.infeas,lastInfeas)>0.95*bestInfeas)
479        majorIterations_ = CoinMin(majorIterations_,iteration); // not getting feasible
480    }
481    lastInfeas = result.infeas;
482    numberAway=n;
483    keepinfeas = result.infeas;
484    lastWeighted=result.weighted;
485    iterationTotal += result.iteration;
486    if (iteration==1) {
487      if ((strategy_&1024)!=0&&mu<1.0e-10) 
488        result.infeas=firstInfeas*0.8;
489      if (majorIterations_>=50)
490        result.infeas *= 0.8;
491      if (result.infeas>firstInfeas*0.9
492          &&result.infeas>reasonableInfeas) {
493        iteration--;
494        if (majorIterations_<50)
495          mu*=1.0e-1;
496        else
497          mu*=0.7;
498        bestFeasible=1.0e31;
499        bestWeighted=1.0e60;
500        numberBaseTrys++;
501        if (mu<1.0e-30||(numberBaseTrys>10&&lightWeight_)) {
502          // back to all slack basis
503          lightWeight_=2;
504          break;
505        }
506        memcpy(colsol,saveSol,ncols*sizeof(double));
507      } else {
508        // Save best solution
509        memcpy(saveSol,colsol,ncols*sizeof(double));
510        maxIts=maxIts2;
511        checkIteration=0;
512        if ((strategy_&1024)!=0) mu *= 1.0e-1;
513      }
514    } else if (result.infeas<bestInfeas) {
515      // Save best solution
516      memcpy(saveSol,colsol,ncols*sizeof(double));
517    }
518    bestInfeas=CoinMin(bestInfeas,result.infeas);
519    if (iteration) {
520      /* this code is in to force it to terminate sometime */
521      double changeMu=factor;
522      if ((saveStrategy&64)!=0) {
523        keepinfeas=0.0; /* switch off ranga's increase */
524        fakeSmall=smallInfeas;
525      } else {
526        fakeSmall=-1.0;
527      }
528      saveLambdaScale=0.0;
529      if (result.infeas>reasonableInfeas||
530          (nTry+1==maxBigIts&&result.infeas>fakeSmall)) {
531        if (result.infeas>lastResult.infeas*(1.0-dropEnoughFeasibility_)||
532            nTry+1==maxBigIts||
533            (result.infeas>lastResult.infeas*0.9
534             &&result.weighted>lastResult.weighted
535             -dropEnoughWeighted_*fabs(lastResult.weighted))) {
536          mu*=changeMu;
537          if ((saveStrategy&32)!=0&&result.infeas<reasonableInfeas&&0) {
538            reasonableInfeas=CoinMax(smallInfeas,reasonableInfeas*sqrt(changeMu));
539            printf("reasonable infeas now %g\n",reasonableInfeas);
540          }
541          nTry=0;
542          bestFeasible=1.0e31;
543          bestWeighted=1.0e60;
544          checkIteration=0;
545          lambdaIteration=0;
546#define LAMBDA
547#ifdef LAMBDA
548          if ((saveStrategy&2048)==0) {
549            memset(lambda,0,nrows*sizeof(double));
550          }
551#else
552          memset(lambda,0,nrows*sizeof(double));
553#endif
554        } else {
555          nTry++;
556        }
557      } else if (lambdaIterations_>=0) {
558        /* update lambda  */
559        double scale=1.0/mu;
560        int i,nnz=0;
561        saveLambdaScale=scale;
562         lambdaIteration++;
563         if ((saveStrategy&4)==0) drop = drop_/50.0;
564         if (lambdaIteration>4 && 
565            (((lambdaIteration%10)==0 && smallInfeas<keepinfeas) ||
566             (lambdaIteration%5)==0 && 1.5*smallInfeas<keepinfeas)) {
567           //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
568           smallInfeas *= 1.5;
569         }
570         if ((saveStrategy&2048)==0) {
571           for (i=0;i<nrows;i++) {
572             if (lambda[i]) nnz++;
573             lambda[i]+= scale*rowsol[i];
574           }
575         } else {
576           nnz=1;
577#ifdef LAMBDA
578           for (i=0;i<nrows;i++) {
579             lambda[i]+= scale*rowsol[i];
580           }
581#else
582           for (i=0;i<nrows;i++) {
583             lambda[i] = scale*rowsol[i];
584           }
585           for (i=0;i<ncols;i++) {
586             CoinBigIndex j;
587             double value=cost[i]*maxmin;
588             for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
589               int irow=row[j];
590               value+=element[j]*lambda[irow];
591             }
592             cost[i]=value*maxmin;
593           }
594           for (i=0;i<nrows;i++) {
595             offset+=lambda[i]*rowupper[i];
596             lambda[i]=0.0;
597           }
598#ifdef DEBUG
599           printf("offset %g\n",offset);
600#endif
601           model_->setDblParam(OsiObjOffset,offset);
602#endif
603         }
604        nTry++;
605        if (!nnz) {
606          bestFeasible=1.0e32;
607          bestWeighted=1.0e60;
608          checkIteration=0;
609          result.weighted=1.0e31;
610        }
611#ifdef DEBUG
612        double trueCost=0.0;
613        for (i=0;i<ncols;i++) {
614          int j;
615          trueCost+=cost[i]*colsol[i];
616        }
617        printf("True objective %g\n",trueCost-offset);
618#endif
619      } else {
620        nTry++;
621      }
622      lastResult=result;
623      if (result.infeas<reasonableInfeas&&!belowReasonable) {
624        belowReasonable=1;
625        bestFeasible=1.0e32;
626        bestWeighted=1.0e60;
627        checkIteration=0;
628        result.weighted=1.0e31;
629      }
630    }
631    if (iteration>=majorIterations_) {
632      // If small and not feasible and crash then dive dive dive
633      if (0&&result.infeas>1.0&&majorIterations_<30&&(maxIts2_==11||maxIts2_==23)) {
634        maxIts=7;
635        majorIterations_=100;
636      } else {
637        if (logLevel>2) 
638          printf("Exiting due to number of major iterations\n");
639        break;
640      }
641    }
642  }
643  majorIterations_ = saveMajorIterations;
644  // put back best solution
645  memcpy(colsol,saveSol,ncols*sizeof(double));
646#ifndef OSI_IDIOT
647  if (scaled) {
648    // Scale solution and free arrays
649    const double * rowScale = model_->rowScale();
650    const double * columnScale = model_->columnScale();
651    int icol,irow;
652    for (icol=0;icol<ncols;icol++) {
653      colsol[icol] *= columnScale[icol];
654      dj[icol] /= columnScale[icol];
655    }
656    for (irow=0;irow<nrows;irow++) {
657      rowsol[irow] /= rowScale[irow];
658      pi[irow] *= rowScale[irow];
659    }
660    // Don't know why getting Microsoft problems
661#if defined (_MSC_VER)
662    delete [] ( double *) elemXX;
663#else
664    delete [] elemXX;
665#endif
666    model_->setRowScale(NULL);
667    model_->setColumnScale(NULL);
668    delete [] lower;
669    delete [] upper;
670    delete [] cost;
671    lower = model_->columnLower();
672    upper = model_->columnUpper();
673    cost = model_->objective();
674    //rowlower = model_->rowLower();
675  }
676#endif
677#define TRYTHIS
678#ifdef TRYTHIS
679  if ((saveStrategy&2048)!=0) {
680    double offset;
681    model_->getDblParam(OsiObjOffset,offset);
682    for (i=0;i<ncols;i++) {
683      CoinBigIndex j;
684      double djval=cost[i]*maxmin;
685      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
686        int irow=row[j];
687        djval -= element[j]*lambda[irow];
688      }
689      cost[i]=djval;
690    }
691    for (i=0;i<nrows;i++) {
692      offset+=lambda[i]*rowupper[i];
693    }
694    model_->setDblParam(OsiObjOffset,offset);
695  }
696#endif
697  if (saveLambdaScale) {
698    /* back off last update */
699    for (i=0;i<nrows;i++) {
700      lambda[i]-= saveLambdaScale*rowsol[i];
701    }
702  }
703  muAtExit_=mu;
704  n=0;
705  for (i=ordStart;i<ordEnd;i++) {
706    if (colsol[i]>lower[i]+fixTolerance) {
707      n++;
708      whenUsed[i]=iteration;
709    } else {
710      colsol[i]=lower[i];
711    }
712  }
713  if ((logLevel&1)==0) {
714    printf(
715            "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
716            iteration,mu,lastResult.infeas,lastResult.objval,n);
717  }
718#ifndef OSI_IDIOT
719  model_->setSumPrimalInfeasibilities(lastResult.infeas);
720#endif
721  {
722    double large=0.0;
723    int i;
724    memset(rowsol,0,nrows*sizeof(double));
725    for (i=0;i<ncols;i++) {
726      CoinBigIndex j;
727      double value=colsol[i];
728      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
729        int irow=row[j];
730        rowsol[irow] += element[j]*value;
731      }
732    }
733    for (i=0;i<nrows;i++) {
734      if (rowsol[i] > rowupper[i]) {
735        double diff=rowsol[i]-rowupper[i];
736        if (diff>large) 
737          large=diff;
738      } else if (rowsol[i] < rowlower[i]) {
739        double diff=rowlower[i]-rowsol[i];
740        if (diff>large) 
741          large=diff;
742      } 
743    }
744    if (logLevel>2)
745      printf("largest infeasibility is %g\n",large);
746  }
747  /* subtract out lambda */
748  for (i=0;i<nrows;i++) {
749    pi[i]-=lambda[i];
750  }
751  for (i=0;i<ncols;i++) {
752    CoinBigIndex j;
753    double djval=cost[i]*maxmin;
754    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
755      int irow=row[j];
756      djval -= element[j]*pi[irow];
757    }
758    dj[i]=djval;
759  }
760  if ((strategy_&1024)!=0) {
761    double ratio = ((double) ncols)/((double) nrows);
762    printf("col/row ratio %g infeas ratio %g\n",ratio,lastResult.infeas/firstInfeas);
763    if (lastResult.infeas>0.01*firstInfeas*ratio) {
764      strategy_ &= (~1024);
765      printf(" - layer off\n");
766    } else {
767      printf(" - layer on\n");
768    }
769  }
770  delete [] saveSol;
771  delete [] lambda;
772  // save solution
773  // duals not much use - but save anyway
774#ifndef OSI_IDIOT
775  memcpy(model_->primalRowSolution(),rowsol,nrows*sizeof(double));
776  memcpy(model_->primalColumnSolution(),colsol,ncols*sizeof(double));
777  memcpy(model_->dualRowSolution(),pi,nrows*sizeof(double));
778  memcpy(model_->dualColumnSolution(),dj,ncols*sizeof(double));
779#else
780  model_->setColSolution(colsol);
781  model_->setRowPrice(pi);
782  delete [] cost;
783#endif
784  delete [] rowsol;
785  delete [] colsol;
786  delete [] pi;
787  delete [] dj;
788  delete [] rowlower;
789  delete [] rowupper;
790  return ;
791}
792#ifndef OSI_IDIOT
793void
794Idiot::crossOver(int mode)
795{
796  if (lightWeight_==2) {
797    // total failure
798    model_->allSlackBasis();
799    return;
800  }
801  double fixTolerance=IDIOT_FIX_TOLERANCE;
802  double startTime = CoinCpuTime();
803  ClpSimplex * saveModel=NULL;
804  ClpMatrixBase * matrix = model_->clpMatrix();
805  const int * row = matrix->getIndices();
806  const CoinBigIndex * columnStart = matrix->getVectorStarts();
807  const int * columnLength = matrix->getVectorLengths(); 
808  const double * element = matrix->getElements();
809  const double * rowupper = model_->getRowUpper();
810  int nrows=model_->getNumRows();
811  int ncols=model_->getNumCols();
812  double * rowsol, * colsol;
813  // different for Osi
814  double * lower = model_->columnLower();
815  double * upper = model_->columnUpper();
816  const double * rowlower= model_->getRowLower();
817  int * whenUsed=whenUsed_;
818  rowsol= model_->primalRowSolution();
819  colsol= model_->primalColumnSolution();;
820  double * cost=model_->objective();
821
822  int slackEnd,ordStart,ordEnd;
823  int slackStart = countCostedSlacks(model_);
824
825  int addAll = mode&7;
826  int presolve=0;
827
828  double djTolerance = djTolerance_;
829  if (djTolerance>0.0&&djTolerance<1.0)
830    djTolerance=1.0;
831  int iteration;
832  int i, n=0;
833  double ratio=1.0;
834  double objValue=0.0;
835  if ((strategy_&128)!=0) {
836    fixTolerance=SMALL_IDIOT_FIX_TOLERANCE;
837  }
838  if ((mode&16)!=0&&addAll<3) presolve=1;
839  double * saveUpper = NULL;
840  double * saveLower = NULL;
841  if (addAll<3) {
842    saveUpper = new double [ncols];
843    saveLower = new double [ncols];
844    memcpy(saveUpper,upper,ncols*sizeof(double));
845    memcpy(saveLower,lower,ncols*sizeof(double));
846  }
847  if (slackStart>=0) {
848    slackEnd=slackStart+nrows;
849    if (slackStart) {
850      ordStart=0;
851      ordEnd=slackStart;
852    } else {
853      ordStart=nrows;
854      ordEnd=ncols;
855    }
856  } else {
857    slackEnd=slackStart;
858    ordStart=0;
859    ordEnd=ncols;
860  }
861  /* get correct rowsol (without known slacks) */
862  memset(rowsol,0,nrows*sizeof(double));
863  for (i=ordStart;i<ordEnd;i++) {
864    CoinBigIndex j;
865    double value=colsol[i];
866    if (value<lower[i]+fixTolerance) {
867      value=lower[i];
868      colsol[i]=value;
869    }
870    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
871      int irow=row[j];
872      rowsol[irow]+=value*element[j];
873    }
874  }
875  if (slackStart>=0) {
876    for (i=0;i<nrows;i++) {
877      if (ratio*rowsol[i]>rowlower[i]&&rowsol[i]>1.0e-8) {
878        ratio=rowlower[i]/rowsol[i];
879      }
880    }
881    for (i=0;i<nrows;i++) {
882      rowsol[i]*=ratio;
883    }
884    for (i=ordStart;i<ordEnd;i++) {
885      double value=colsol[i]*ratio;
886      colsol[i]=value;
887      objValue+=value*cost[i];
888    }
889    for (i=0;i<nrows;i++) {
890      double value=rowlower[i]-rowsol[i];
891      colsol[i+slackStart]=value;
892      objValue+=value*cost[i+slackStart];
893    }
894    printf("New objective after scaling %g\n",objValue);
895  }
896#if 0
897   maybe put back - but just get feasible ?
898  // If not many fixed then just exit
899  int numberFixed=0;
900  for (i=ordStart;i<ordEnd;i++) {
901    if (colsol[i]<lower[i]+fixTolerance)
902      numberFixed++;
903    else if (colsol[i]>upper[i]-fixTolerance)
904      numberFixed++;
905  }
906  if (numberFixed<ncols/2) {
907    addAll=3;
908    presolve=0;
909  }
910#endif
911  model_->createStatus();
912  /* addAll
913     0 - chosen,all used, all
914     1 - chosen, all
915     2 - all
916     3 - do not do anything  - maybe basis
917  */
918  for (i=ordStart;i<ordEnd;i++) {
919    if (addAll<2) {
920      if (colsol[i]<lower[i]+fixTolerance) {
921        upper[i]=lower[i];
922        colsol[i]=lower[i];
923      } else if (colsol[i]>upper[i]-fixTolerance) {
924        lower[i]=upper[i];
925        colsol[i]=upper[i];
926      }
927    }
928    model_->setColumnStatus(i,ClpSimplex::superBasic);
929  }
930  double maxmin;
931  if (model_->getObjSense()==-1.0) {
932    maxmin=-1.0;
933  } else {
934    maxmin=1.0;
935  }
936  if (slackStart>=0) {
937    for (i=0;i<nrows;i++) {
938      model_->setRowStatus(i,ClpSimplex::superBasic);
939    }
940    for (i=slackStart;i<slackEnd;i++) {
941      model_->setColumnStatus(i,ClpSimplex::basic);
942    }
943  } else {
944    /* still try and put singletons rather than artificials in basis */
945    int ninbas=0;
946    for (i=0;i<nrows;i++) {
947      model_->setRowStatus(i,ClpSimplex::basic);
948    }
949    for (i=0;i<ncols;i++) {
950      if (columnLength[i]==1&&upper[i]>lower[i]+1.0e-5) {
951        CoinBigIndex j =columnStart[i];
952        double value=element[j];
953        int irow=row[j];
954        double rlo=rowlower[irow];
955        double rup=rowupper[irow];
956        double clo=lower[i];
957        double cup=upper[i];
958        double csol=colsol[i];
959        /* adjust towards feasibility */
960        double move=0.0;
961        if (rowsol[irow]>rup) {
962          move=(rup-rowsol[irow])/value;
963          if (value>0.0) {
964            /* reduce */
965            if (csol+move<clo) move=CoinMin(0.0,clo-csol);
966          } else {
967            /* increase */
968            if (csol+move>cup) move=CoinMax(0.0,cup-csol);
969          }
970        } else if (rowsol[irow]<rlo) {
971          move=(rlo-rowsol[irow])/value;
972          if (value>0.0) {
973            /* increase */
974            if (csol+move>cup) move=CoinMax(0.0,cup-csol);
975          } else {
976            /* reduce */
977            if (csol+move<clo) move=CoinMin(0.0,clo-csol);
978          }
979        } else {
980          /* move to improve objective */
981          if (cost[i]*maxmin>0.0) {
982            if (value>0.0) {
983              move=(rlo-rowsol[irow])/value;
984              /* reduce */
985              if (csol+move<clo) move=CoinMin(0.0,clo-csol);
986            } else {
987              move=(rup-rowsol[irow])/value;
988              /* increase */
989              if (csol+move>cup) move=CoinMax(0.0,cup-csol);
990            }
991          } else if (cost[i]*maxmin<0.0) {
992            if (value>0.0) {
993              move=(rup-rowsol[irow])/value;
994              /* increase */
995              if (csol+move>cup) move=CoinMax(0.0,cup-csol);
996            } else {
997              move=(rlo-rowsol[irow])/value;
998              /* reduce */
999              if (csol+move<clo) move=CoinMin(0.0,clo-csol);
1000            }
1001          }
1002        }
1003        rowsol[irow] +=move*value;
1004        colsol[i]+=move;
1005        /* put in basis if row was artificial */
1006        if (rup-rlo<1.0e-7&&model_->getRowStatus(irow)==ClpSimplex::basic) {
1007          model_->setRowStatus(irow,ClpSimplex::superBasic);
1008          model_->setColumnStatus(i,ClpSimplex::basic);
1009          ninbas++;
1010        }
1011      }
1012    }
1013    /*printf("%d in basis\n",ninbas);*/
1014  }
1015  if (addAll<3) {
1016    ClpPresolve pinfo;
1017    if (presolve) {
1018      saveModel = model_;
1019      model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1020    }
1021    if (model_) {
1022      model_->primal(1);
1023      if (presolve) {
1024        pinfo.postsolve(true);
1025        delete model_;
1026        model_ = saveModel;
1027        saveModel=NULL;
1028      }
1029    } else {
1030      // not feasible
1031      addAll=1;
1032      presolve=0;
1033      model_ = saveModel;
1034      saveModel=NULL;
1035    }
1036    if (addAll<2) {
1037      n=0;
1038      if (!addAll ) {
1039        /* could do scans to get a good number */
1040        iteration=1;
1041        for (i=ordStart;i<ordEnd;i++) {
1042          if (whenUsed[i]>=iteration) {
1043            if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1044              n++;
1045              upper[i]=saveUpper[i];
1046              lower[i]=saveLower[i];
1047            }
1048          }
1049        }
1050      } else {
1051        for (i=ordStart;i<ordEnd;i++) {
1052          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1053            n++;
1054            upper[i]=saveUpper[i];
1055            lower[i]=saveLower[i];
1056          }
1057        }
1058      }
1059      printf("Time so far %g, %d now added from previous iterations\n",
1060             CoinCpuTime()-startTime,n);
1061      if (addAll)
1062        presolve=0;
1063      if (presolve) {
1064        saveModel = model_;
1065        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1066      } else {
1067        presolve=0;
1068      }
1069      model_->primal(1);
1070      if (presolve) {
1071        pinfo.postsolve(true);
1072        delete model_;
1073        model_ = saveModel;
1074        saveModel=NULL;
1075      }
1076      if (!addAll) {
1077        n=0;
1078        for (i=ordStart;i<ordEnd;i++) {
1079          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1080            n++;
1081            upper[i]=saveUpper[i];
1082            lower[i]=saveLower[i];
1083          }
1084        }
1085        printf("Time so far %g, %d now added from previous iterations\n",
1086               CoinCpuTime()-startTime,n);
1087      }
1088      if (presolve) {
1089        saveModel = model_;
1090        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1091      } else {
1092        presolve=0;
1093      }
1094      model_->primal(1);
1095      if (presolve) {
1096        pinfo.postsolve(true);
1097        delete model_;
1098        model_ = saveModel;
1099        saveModel=NULL;
1100      }
1101    }
1102    printf("Total time in crossover %g\n", CoinCpuTime()-startTime);
1103    delete [] saveUpper;
1104    delete [] saveLower;
1105  }
1106  return ;
1107}
1108#endif
1109/*****************************************************************************/
1110
1111// Default contructor
1112Idiot::Idiot()
1113{
1114  model_ = NULL;
1115  maxBigIts_ = 3;
1116  maxIts_ = 5;
1117  logLevel_ = 1; 
1118  logFreq_ = 100;
1119  maxIts2_ = 100;
1120  djTolerance_ = 1e-1;
1121  mu_ = 1e-4;
1122  drop_ = 5.0;
1123  exitDrop_=-1.0e20;
1124  muFactor_ = 0.3333;
1125  stopMu_ = 1e-12;
1126  smallInfeas_ = 1e-1;
1127  reasonableInfeas_ = 1e2;
1128  muAtExit_ =1.0e31;
1129  strategy_ =8;
1130  lambdaIterations_ =0;
1131  checkFrequency_ =100;
1132  whenUsed_ = NULL;
1133  majorIterations_ =30;
1134  exitFeasibility_ =-1.0;
1135  dropEnoughFeasibility_ =0.02;
1136  dropEnoughWeighted_ =0.01;
1137  // adjust
1138  double nrows=10000.0;
1139  int baseIts =(int) sqrt((double)nrows);
1140  baseIts =baseIts/10;
1141  baseIts *= 10;
1142  maxIts2_ =200+baseIts+5;
1143  reasonableInfeas_ =((double) nrows)*0.05;
1144  lightWeight_=0;
1145}
1146// Constructor from model
1147Idiot::Idiot(OsiSolverInterface &model)
1148{
1149  model_ = & model;
1150  maxBigIts_ = 3;
1151  maxIts_ = 5;
1152  logLevel_ = 1; 
1153  logFreq_ = 100;
1154  maxIts2_ = 100;
1155  djTolerance_ = 1e-1;
1156  mu_ = 1e-4;
1157  drop_ = 5.0;
1158  exitDrop_=-1.0e20;
1159  muFactor_ = 0.3333;
1160  stopMu_ = 1e-12;
1161  smallInfeas_ = 1e-1;
1162  reasonableInfeas_ = 1e2;
1163  muAtExit_ =1.0e31;
1164  strategy_ =8;
1165  lambdaIterations_ =0;
1166  checkFrequency_ =100;
1167  whenUsed_ = NULL;
1168  majorIterations_ =30;
1169  exitFeasibility_ =-1.0;
1170  dropEnoughFeasibility_ =0.02;
1171  dropEnoughWeighted_ =0.01;
1172  // adjust
1173  double nrows;
1174  if (model_)
1175    nrows=model_->getNumRows();
1176  else
1177    nrows=10000.0;
1178  int baseIts =(int) sqrt((double)nrows);
1179  baseIts =baseIts/10;
1180  baseIts *= 10;
1181  maxIts2_ =200+baseIts+5;
1182  reasonableInfeas_ =((double) nrows)*0.05;
1183  lightWeight_=0;
1184}
1185// Copy constructor.
1186Idiot::Idiot(const Idiot &rhs)
1187{
1188  model_ = rhs.model_;
1189  if (model_&&rhs.whenUsed_) {
1190    int numberColumns = model_->getNumCols();
1191    whenUsed_ = new int [numberColumns];
1192    memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1193  } else {
1194    whenUsed_=NULL;
1195  }
1196  djTolerance_ = rhs.djTolerance_;
1197  mu_ = rhs.mu_;
1198  drop_ = rhs.drop_;
1199  muFactor_ = rhs.muFactor_;
1200  stopMu_ = rhs.stopMu_;
1201  smallInfeas_ = rhs.smallInfeas_;
1202  reasonableInfeas_ = rhs.reasonableInfeas_;
1203  exitDrop_ = rhs.exitDrop_;
1204  muAtExit_ = rhs.muAtExit_;
1205  exitFeasibility_ = rhs.exitFeasibility_;
1206  dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1207  dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1208  maxBigIts_ = rhs.maxBigIts_;
1209  maxIts_ = rhs.maxIts_;
1210  majorIterations_ = rhs.majorIterations_;
1211  logLevel_ = rhs.logLevel_;
1212  logFreq_ = rhs.logFreq_;
1213  checkFrequency_ = rhs.checkFrequency_;
1214  lambdaIterations_ = rhs.lambdaIterations_;
1215  maxIts2_ = rhs.maxIts2_;
1216  strategy_ = rhs.strategy_;
1217  lightWeight_=rhs.lightWeight_;
1218}
1219// Assignment operator. This copies the data
1220Idiot & 
1221Idiot::operator=(const Idiot & rhs)
1222{
1223  if (this != &rhs) {
1224    delete [] whenUsed_;
1225    model_ = rhs.model_;
1226    if (model_&&rhs.whenUsed_) {
1227      int numberColumns = model_->getNumCols();
1228      whenUsed_ = new int [numberColumns];
1229      memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1230    } else {
1231      whenUsed_=NULL;
1232    }
1233    djTolerance_ = rhs.djTolerance_;
1234    mu_ = rhs.mu_;
1235    drop_ = rhs.drop_;
1236    muFactor_ = rhs.muFactor_;
1237    stopMu_ = rhs.stopMu_;
1238    smallInfeas_ = rhs.smallInfeas_;
1239    reasonableInfeas_ = rhs.reasonableInfeas_;
1240    exitDrop_ = rhs.exitDrop_;
1241    muAtExit_ = rhs.muAtExit_;
1242    exitFeasibility_ = rhs.exitFeasibility_;
1243    dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1244    dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1245    maxBigIts_ = rhs.maxBigIts_;
1246    maxIts_ = rhs.maxIts_;
1247    majorIterations_ = rhs.majorIterations_;
1248    logLevel_ = rhs.logLevel_;
1249    logFreq_ = rhs.logFreq_;
1250    checkFrequency_ = rhs.checkFrequency_;
1251    lambdaIterations_ = rhs.lambdaIterations_;
1252    maxIts2_ = rhs.maxIts2_;
1253    strategy_ = rhs.strategy_;
1254    lightWeight_=rhs.lightWeight_;
1255  }
1256  return *this;
1257}
1258Idiot::~Idiot()
1259{
1260  delete [] whenUsed_;
1261}
Note: See TracBrowser for help on using the repository browser.