source: trunk/Idiot.cpp @ 399

Last change on this file since 399 was 399, checked in by forrest, 16 years ago

to CoinMax? and CoinMin?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 32.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 "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  double * rowlower= model_->rowLower();
228#else
229  double * cost = new double [ncols];
230  memcpy( cost, model_->getObjCoefficients(), ncols*sizeof(double));
231  const double * lower = model_->getColLower();
232  const double * upper = model_->getColUpper();
233  const double * rowlower= model_->getRowLower();
234#endif
235  const double *elemXX;
236  double * saveSol;
237  double * rowupper= new double[nrows]; // not const as modified
238  memcpy(rowupper,model_->getRowUpper(),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      rowlower = new double[nrows];
326      memcpy(rowlower,model_->rowLower(),nrows*sizeof(double));
327      for (irow=0;irow<nrows;irow++) {
328        double multiplier = rowScale[irow];
329        if (rowlower[irow]>-1.0e50)
330          rowlower[irow] *= multiplier;
331        if (rowupper[irow]<1.0e50)
332          rowupper[irow] *= multiplier;
333        rowsol[irow] *= multiplier;
334      }
335      int length = columnStart[ncols-1]+columnLength[ncols-1];
336      double * elemYY = new double[length];
337      for (i=0;i<ncols;i++) {
338        CoinBigIndex j;
339        double scale = columnScale[i];
340        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
341          int irow=row[j];
342          elemYY[j] = element[j]*scale*rowScale[irow];
343        }
344      }
345      elemXX=elemYY;
346    }
347  }
348#endif
349  for (i=0;i<ncols;i++) {
350    CoinBigIndex j;
351    double dd=columnLength[i];
352    dd=cost[i]/dd;
353    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
354      int irow=row[j];
355      if (dd<rowsol[irow]) {
356        rowsol[irow]=dd;
357      }
358    }
359  }
360  d2=0.0;
361  for (i=0;i<nrows;i++) {
362    d2+=rowsol[i];
363  }
364  d2*=2.0; /* for luck */
365 
366  d2=d2/((double) (4*nrows+8000));
367  d2*=0.5; /* halve with more flexible method */
368  if (d2<5.0) d2=5.0;
369  if (djExit==0.0) {
370    djExit=d2;
371  }
372  if ((saveStrategy&4)!=0) {
373    /* go to relative tolerances - first small */
374    djExit=1.0e-10;
375    djFlag=1.0e-5;
376    drop=1.0e-10;
377  }
378  memset(whenUsed,0,ncols*sizeof(int));
379  strategy=strategies[strategy];
380  if ((saveStrategy&8)!=0) strategy |= 64; /* don't allow large theta */
381  memcpy(saveSol,colsol,ncols*sizeof(double));
382 
383  lastResult=IdiSolve(nrows,ncols,rowsol ,colsol,pi,
384                       dj,cost,rowlower,rowupper,
385                       lower,upper,elemXX,row,columnStart,columnLength,lambda,
386                       0,mu,drop,
387                       maxmin,offset,strategy,djTol,djExit,djFlag);
388  n=0;
389  for (i=ordStart;i<ordEnd;i++) {
390    if (colsol[i]>lower[i]+fixTolerance) {
391      if (colsol[i]<upper[i]-fixTolerance) {
392        n++;
393      } else {
394        colsol[i]=upper[i];
395      }
396      whenUsed[i]=iteration;
397    } else {
398      colsol[i]=lower[i];
399    }
400  }
401  if ((logLevel_&1)!=0) {
402#ifndef OSI_IDIOT
403    if (!handler) {
404#endif
405      printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
406             iteration,lastResult.infeas,lastResult.objval,mu,lastResult.iteration,n);
407#ifndef OSI_IDIOT
408    } else {
409      handler->message(CLP_IDIOT_ITERATION,*messages)
410        <<iteration<<lastResult.infeas<<lastResult.objval<<mu<<lastResult.iteration<<n
411        <<CoinMessageEol;
412    }
413#endif
414  }
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    bestInfeas=CoinMin(bestInfeas,result.infeas);
482    lastInfeas = result.infeas;
483    numberAway=n;
484    keepinfeas = result.infeas;
485    lastWeighted=result.weighted;
486    iterationTotal += result.iteration;
487    if (iteration==1) {
488      if ((strategy_&1024)!=0&&mu<1.0e-10) 
489        result.infeas=firstInfeas*0.8;
490      if (result.infeas>firstInfeas*0.9
491          &&result.infeas>reasonableInfeas) {
492        iteration--;
493        memcpy(colsol,saveSol,ncols*sizeof(double));
494        if (majorIterations_<50)
495          mu*=1.0e-1;
496        else
497          mu*=0.7;
498        bestFeasible=1.0e31;
499        bestWeighted=1.0e60;
500        if (mu<1.0e-30) break;
501      } else {
502        delete [] saveSol;
503        saveSol=0;
504        maxIts=maxIts2;
505        checkIteration=0;
506        if ((strategy_&1024)!=0) mu *= 1.0e-1;
507      }
508    }
509    if (iteration) {
510      /* this code is in to force it to terminate sometime */
511      double changeMu=factor;
512      if ((saveStrategy&64)!=0) {
513        keepinfeas=0.0; /* switch off ranga's increase */
514        fakeSmall=smallInfeas;
515      } else {
516        fakeSmall=-1.0;
517      }
518      saveLambdaScale=0.0;
519      if (result.infeas>reasonableInfeas||
520          (nTry+1==maxBigIts&&result.infeas>fakeSmall)) {
521        if (result.infeas>lastResult.infeas*(1.0-dropEnoughFeasibility_)||
522            nTry+1==maxBigIts||
523            (result.infeas>lastResult.infeas*0.9
524             &&result.weighted>lastResult.weighted
525             -dropEnoughWeighted_*fabs(lastResult.weighted))) {
526          mu*=changeMu;
527          if ((saveStrategy&32)!=0&&result.infeas<reasonableInfeas&&0) {
528            reasonableInfeas=CoinMax(smallInfeas,reasonableInfeas*sqrt(changeMu));
529            printf("reasonable infeas now %g\n",reasonableInfeas);
530          }
531          nTry=0;
532          bestFeasible=1.0e31;
533          bestWeighted=1.0e60;
534          checkIteration=0;
535          lambdaIteration=0;
536#define LAMBDA
537#ifdef LAMBDA
538          if ((saveStrategy&2048)==0) {
539            memset(lambda,0,nrows*sizeof(double));
540          }
541#else
542          memset(lambda,0,nrows*sizeof(double));
543#endif
544        } else {
545          nTry++;
546        }
547      } else if (lambdaIterations_>=0) {
548        /* update lambda  */
549        double scale=1.0/mu;
550        int i,nnz=0;
551        saveLambdaScale=scale;
552         lambdaIteration++;
553         if ((saveStrategy&4)==0) drop = drop_/50.0;
554         if (lambdaIteration>4 && 
555            (((lambdaIteration%10)==0 && smallInfeas<keepinfeas) ||
556             (lambdaIteration%5)==0 && 1.5*smallInfeas<keepinfeas)) {
557           //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
558           smallInfeas *= 1.5;
559         }
560         if ((saveStrategy&2048)==0) {
561           for (i=0;i<nrows;i++) {
562             if (lambda[i]) nnz++;
563             lambda[i]+= scale*rowsol[i];
564           }
565         } else {
566           nnz=1;
567#ifdef LAMBDA
568           for (i=0;i<nrows;i++) {
569             lambda[i]+= scale*rowsol[i];
570           }
571#else
572           for (i=0;i<nrows;i++) {
573             lambda[i] = scale*rowsol[i];
574           }
575           for (i=0;i<ncols;i++) {
576             CoinBigIndex j;
577             double value=cost[i]*maxmin;
578             for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
579               int irow=row[j];
580               value+=element[j]*lambda[irow];
581             }
582             cost[i]=value*maxmin;
583           }
584           for (i=0;i<nrows;i++) {
585             offset+=lambda[i]*rowupper[i];
586             lambda[i]=0.0;
587           }
588#ifdef DEBUG
589           printf("offset %g\n",offset);
590#endif
591           model_->setDblParam(OsiObjOffset,offset);
592#endif
593         }
594        nTry++;
595        if (!nnz) {
596          bestFeasible=1.0e32;
597          bestWeighted=1.0e60;
598          checkIteration=0;
599          result.weighted=1.0e31;
600        }
601#ifdef DEBUG
602        for (i=0;i<ncols;i++) {
603          int j;
604          trueCost+=cost[i]*colsol[i];
605        }
606        printf("True objective %g\n",trueCost-offset);
607#endif
608      } else {
609        nTry++;
610      }
611      lastResult=result;
612      if (result.infeas<reasonableInfeas&&!belowReasonable) {
613        belowReasonable=1;
614        bestFeasible=1.0e32;
615        bestWeighted=1.0e60;
616        checkIteration=0;
617        result.weighted=1.0e31;
618      }
619    }
620    if (iteration>=majorIterations_) {
621      // If small and not feasible and crash then dive dive dive
622      if (0&&result.infeas>1.0&&majorIterations_<30&&(maxIts2_==11||maxIts2_==23)) {
623        maxIts=7;
624        majorIterations_=100;
625      } else {
626        if (logLevel>2) 
627          printf("Exiting due to number of major iterations\n");
628        break;
629      }
630    }
631  }
632  majorIterations_ = saveMajorIterations;
633#ifndef OSI_IDIOT
634  if (scaled) {
635    // Scale solution and free arrays
636    const double * rowScale = model_->rowScale();
637    const double * columnScale = model_->columnScale();
638    int icol,irow;
639    for (icol=0;icol<ncols;icol++) {
640      colsol[icol] *= columnScale[icol];
641      dj[icol] /= columnScale[icol];
642    }
643    for (irow=0;irow<nrows;irow++) {
644      rowsol[irow] /= rowScale[irow];
645      pi[irow] *= rowScale[irow];
646    }
647    // Don't know why getting Microsoft problems
648#if defined (_MSC_VER)
649    delete [] ( double *)rowScale;
650    delete [] ( double *) columnScale;
651    delete [] ( double *) elemXX;
652#else
653    delete [] rowScale;
654    delete [] columnScale;
655    delete [] elemXX;
656#endif
657    model_->setRowScale(NULL);
658    model_->setColumnScale(NULL);
659    delete [] lower;
660    delete [] upper;
661    delete [] cost;
662    lower = model_->columnLower();
663    upper = model_->columnUpper();
664    cost = model_->objective();
665    rowlower = model_->rowLower();
666  }
667#endif
668#define TRYTHIS
669#ifdef TRYTHIS
670  if ((saveStrategy&2048)!=0) {
671    double offset;
672    model_->getDblParam(OsiObjOffset,offset);
673    for (i=0;i<ncols;i++) {
674      CoinBigIndex j;
675      double djval=cost[i]*maxmin;
676      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
677        int irow=row[j];
678        djval -= element[j]*lambda[irow];
679      }
680      cost[i]=djval;
681    }
682    for (i=0;i<nrows;i++) {
683      offset+=lambda[i]*rowupper[i];
684    }
685    model_->setDblParam(OsiObjOffset,offset);
686  }
687#endif
688  if (saveLambdaScale) {
689    /* back off last update */
690    for (i=0;i<nrows;i++) {
691      lambda[i]-= saveLambdaScale*rowsol[i];
692    }
693  }
694  muAtExit_=mu;
695  n=0;
696  for (i=ordStart;i<ordEnd;i++) {
697    if (colsol[i]>lower[i]+fixTolerance) {
698      n++;
699      whenUsed[i]=iteration;
700    } else {
701      colsol[i]=lower[i];
702    }
703  }
704  if ((logLevel&1)==0) {
705    printf(
706            "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
707            iteration,mu,lastResult.infeas,lastResult.objval,n);
708  }
709#ifndef OSI_IDIOT
710  model_->setSumPrimalInfeasibilities(lastResult.infeas);
711#endif
712  {
713    double large=0.0;
714    int i;
715    memset(rowsol,0,nrows*sizeof(double));
716    for (i=0;i<ncols;i++) {
717      CoinBigIndex j;
718      double value=colsol[i];
719      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
720        int irow=row[j];
721        rowsol[irow] += element[j]*value;
722      }
723    }
724    for (i=0;i<nrows;i++) {
725      if (rowsol[i] > rowupper[i]) {
726        double diff=rowsol[i]-rowupper[i];
727        if (diff>large) 
728          large=diff;
729      } else if (rowsol[i] < rowlower[i]) {
730        double diff=rowlower[i]-rowsol[i];
731        if (diff>large) 
732          large=diff;
733      } 
734    }
735    if (logLevel>2)
736      printf("largest infeasibility is %g\n",large);
737  }
738  /* subtract out lambda */
739  for (i=0;i<nrows;i++) {
740    pi[i]-=lambda[i];
741  }
742  for (i=0;i<ncols;i++) {
743    CoinBigIndex j;
744    double djval=cost[i]*maxmin;
745    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
746      int irow=row[j];
747      djval -= element[j]*pi[irow];
748    }
749    dj[i]=djval;
750  }
751  if ((strategy_&1024)!=0) {
752    double ratio = ((double) ncols)/((double) nrows);
753    printf("col/row ratio %g infeas ratio %g\n",ratio,lastResult.infeas/firstInfeas);
754    if (lastResult.infeas>0.01*firstInfeas*ratio) {
755      strategy_ &= (~1024);
756      printf(" - layer off\n");
757    } else {
758      printf(" - layer on\n");
759    }
760  }
761  delete [] saveSol;
762  delete [] lambda;
763  // save solution
764  // duals not much use - but save anyway
765#ifndef OSI_IDIOT
766  memcpy(model_->primalRowSolution(),rowsol,nrows*sizeof(double));
767  memcpy(model_->primalColumnSolution(),colsol,ncols*sizeof(double));
768  memcpy(model_->dualRowSolution(),pi,nrows*sizeof(double));
769  memcpy(model_->dualColumnSolution(),dj,ncols*sizeof(double));
770#else
771  model_->setColSolution(colsol);
772  model_->setRowPrice(pi);
773  delete [] cost;
774#endif
775  delete [] rowsol;
776  delete [] colsol;
777  delete [] pi;
778  delete [] dj;
779  return ;
780}
781#ifndef OSI_IDIOT
782void
783Idiot::crossOver(int mode)
784{
785  double fixTolerance=IDIOT_FIX_TOLERANCE;
786  double startTime = CoinCpuTime();
787  ClpSimplex * saveModel=NULL;
788  ClpMatrixBase * matrix = model_->clpMatrix();
789  const int * row = matrix->getIndices();
790  const CoinBigIndex * columnStart = matrix->getVectorStarts();
791  const int * columnLength = matrix->getVectorLengths(); 
792  const double * element = matrix->getElements();
793  const double * rowupper = model_->getRowUpper();
794  int nrows=model_->getNumRows();
795  int ncols=model_->getNumCols();
796  double * rowsol, * colsol;
797  // different for Osi
798  double * lower = model_->columnLower();
799  double * upper = model_->columnUpper();
800  const double * rowlower= model_->getRowLower();
801  int * whenUsed=whenUsed_;
802  rowsol= model_->primalRowSolution();
803  colsol= model_->primalColumnSolution();;
804  double * cost=model_->objective();
805
806  int slackEnd,ordStart,ordEnd;
807  int slackStart = countCostedSlacks(model_);
808
809  int addAll = mode&7;
810  int presolve=0;
811
812  double djTolerance = djTolerance_;
813  if (djTolerance>0.0&&djTolerance<1.0)
814    djTolerance=1.0;
815  int iteration;
816  int i, n=0;
817  double ratio=1.0;
818  double objValue=0.0;
819  if ((strategy_&128)!=0) {
820    fixTolerance=SMALL_IDIOT_FIX_TOLERANCE;
821  }
822  if ((mode&16)!=0&&addAll<3) presolve=1;
823  double * saveUpper = NULL;
824  double * saveLower = NULL;
825  if (addAll<3) {
826    saveUpper = new double [ncols];
827    saveLower = new double [ncols];
828    memcpy(saveUpper,upper,ncols*sizeof(double));
829    memcpy(saveLower,lower,ncols*sizeof(double));
830  }
831  if (slackStart>=0) {
832    slackEnd=slackStart+nrows;
833    if (slackStart) {
834      ordStart=0;
835      ordEnd=slackStart;
836    } else {
837      ordStart=nrows;
838      ordEnd=ncols;
839    }
840  } else {
841    slackEnd=slackStart;
842    ordStart=0;
843    ordEnd=ncols;
844  }
845  /* get correct rowsol (without known slacks) */
846  memset(rowsol,0,nrows*sizeof(double));
847  for (i=ordStart;i<ordEnd;i++) {
848    CoinBigIndex j;
849    double value=colsol[i];
850    if (value<lower[i]+fixTolerance) {
851      value=lower[i];
852      colsol[i]=value;
853    }
854    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
855      int irow=row[j];
856      rowsol[irow]+=value*element[j];
857    }
858  }
859  if (slackStart>=0) {
860    for (i=0;i<nrows;i++) {
861      if (ratio*rowsol[i]>rowlower[i]&&rowsol[i]>1.0e-8) {
862        ratio=rowlower[i]/rowsol[i];
863      }
864    }
865    for (i=0;i<nrows;i++) {
866      rowsol[i]*=ratio;
867    }
868    for (i=ordStart;i<ordEnd;i++) {
869      double value=colsol[i]*ratio;
870      colsol[i]=value;
871      objValue+=value*cost[i];
872    }
873    for (i=0;i<nrows;i++) {
874      double value=rowlower[i]-rowsol[i];
875      colsol[i+slackStart]=value;
876      objValue+=value*cost[i+slackStart];
877    }
878    printf("New objective after scaling %g\n",objValue);
879  }
880#if 0
881   maybe put back - but just get feasible ?
882  // If not many fixed then just exit
883  int numberFixed=0;
884  for (i=ordStart;i<ordEnd;i++) {
885    if (colsol[i]<lower[i]+fixTolerance)
886      numberFixed++;
887    else if (colsol[i]>upper[i]-fixTolerance)
888      numberFixed++;
889  }
890  if (numberFixed<ncols/2) {
891    addAll=3;
892    presolve=0;
893  }
894#endif
895  model_->createStatus();
896  /* addAll
897     0 - chosen,all used, all
898     1 - chosen, all
899     2 - all
900     3 - do not do anything  - maybe basis
901  */
902  for (i=ordStart;i<ordEnd;i++) {
903    if (addAll<2) {
904      if (colsol[i]<lower[i]+fixTolerance) {
905        upper[i]=lower[i];
906        colsol[i]=lower[i];
907      } else if (colsol[i]>upper[i]-fixTolerance) {
908        lower[i]=upper[i];
909        colsol[i]=upper[i];
910      }
911    }
912    model_->setColumnStatus(i,ClpSimplex::superBasic);
913  }
914  double maxmin;
915  if (model_->getObjSense()==-1.0) {
916    maxmin=-1.0;
917  } else {
918    maxmin=1.0;
919  }
920  if (slackStart>=0) {
921    for (i=0;i<nrows;i++) {
922      model_->setRowStatus(i,ClpSimplex::superBasic);
923    }
924    for (i=slackStart;i<slackEnd;i++) {
925      model_->setColumnStatus(i,ClpSimplex::basic);
926    }
927  } else {
928    /* still try and put singletons rather than artificials in basis */
929    int ninbas=0;
930    for (i=0;i<nrows;i++) {
931      model_->setRowStatus(i,ClpSimplex::basic);
932    }
933    for (i=0;i<ncols;i++) {
934      if (columnLength[i]==1&&upper[i]>lower[i]+1.0e-5) {
935        CoinBigIndex j =columnStart[i];
936        double value=element[j];
937        int irow=row[j];
938        double rlo=rowlower[irow];
939        double rup=rowupper[irow];
940        double clo=lower[i];
941        double cup=upper[i];
942        double csol=colsol[i];
943        /* adjust towards feasibility */
944        double move=0.0;
945        if (rowsol[irow]>rup) {
946          move=(rup-rowsol[irow])/value;
947          if (value>0.0) {
948            /* reduce */
949            if (csol+move<clo) move=CoinMin(0.0,clo-csol);
950          } else {
951            /* increase */
952            if (csol+move>cup) move=CoinMax(0.0,cup-csol);
953          }
954        } else if (rowsol[irow]<rlo) {
955          move=(rlo-rowsol[irow])/value;
956          if (value>0.0) {
957            /* increase */
958            if (csol+move>cup) move=CoinMax(0.0,cup-csol);
959          } else {
960            /* reduce */
961            if (csol+move<clo) move=CoinMin(0.0,clo-csol);
962          }
963        } else {
964          /* move to improve objective */
965          if (cost[i]*maxmin>0.0) {
966            if (value>0.0) {
967              move=(rlo-rowsol[irow])/value;
968              /* reduce */
969              if (csol+move<clo) move=CoinMin(0.0,clo-csol);
970            } else {
971              move=(rup-rowsol[irow])/value;
972              /* increase */
973              if (csol+move>cup) move=CoinMax(0.0,cup-csol);
974            }
975          } else if (cost[i]*maxmin<0.0) {
976            if (value>0.0) {
977              move=(rup-rowsol[irow])/value;
978              /* increase */
979              if (csol+move>cup) move=CoinMax(0.0,cup-csol);
980            } else {
981              move=(rlo-rowsol[irow])/value;
982              /* reduce */
983              if (csol+move<clo) move=CoinMin(0.0,clo-csol);
984            }
985          }
986        }
987        rowsol[irow] +=move*value;
988        colsol[i]+=move;
989        /* put in basis if row was artificial */
990        if (rup-rlo<1.0e-7&&model_->getRowStatus(irow)==ClpSimplex::basic) {
991          model_->setRowStatus(irow,ClpSimplex::superBasic);
992          model_->setColumnStatus(i,ClpSimplex::basic);
993          ninbas++;
994        }
995      }
996    }
997    /*printf("%d in basis\n",ninbas);*/
998  }
999  if (addAll<3) {
1000    ClpPresolve pinfo;
1001    if (presolve) {
1002      saveModel = model_;
1003      model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1004    }
1005    if (model_) {
1006      model_->primal(1);
1007      if (presolve) {
1008        pinfo.postsolve(true);
1009        delete model_;
1010        model_ = saveModel;
1011        saveModel=NULL;
1012      }
1013    } else {
1014      // not feasible
1015      addAll=1;
1016      presolve=0;
1017      model_ = saveModel;
1018      saveModel=NULL;
1019    }
1020    if (addAll<2) {
1021      n=0;
1022      if (!addAll ) {
1023        /* could do scans to get a good number */
1024        iteration=1;
1025        for (i=ordStart;i<ordEnd;i++) {
1026          if (whenUsed[i]>=iteration) {
1027            if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1028              n++;
1029              upper[i]=saveUpper[i];
1030              lower[i]=saveLower[i];
1031            }
1032          }
1033        }
1034      } else {
1035        for (i=ordStart;i<ordEnd;i++) {
1036          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1037            n++;
1038            upper[i]=saveUpper[i];
1039            lower[i]=saveLower[i];
1040          }
1041        }
1042      }
1043      printf("Time so far %g, %d now added from previous iterations\n",
1044             CoinCpuTime()-startTime,n);
1045      if (addAll)
1046        presolve=0;
1047      if (presolve) {
1048        saveModel = model_;
1049        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1050      } else {
1051        presolve=0;
1052      }
1053      model_->primal(1);
1054      if (presolve) {
1055        pinfo.postsolve(true);
1056        delete model_;
1057        model_ = saveModel;
1058        saveModel=NULL;
1059      }
1060      if (!addAll) {
1061        n=0;
1062        for (i=ordStart;i<ordEnd;i++) {
1063          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
1064            n++;
1065            upper[i]=saveUpper[i];
1066            lower[i]=saveLower[i];
1067          }
1068        }
1069        printf("Time so far %g, %d now added from previous iterations\n",
1070               CoinCpuTime()-startTime,n);
1071      }
1072      if (presolve) {
1073        saveModel = model_;
1074        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
1075      } else {
1076        presolve=0;
1077      }
1078      model_->primal(1);
1079      if (presolve) {
1080        pinfo.postsolve(true);
1081        delete model_;
1082        model_ = saveModel;
1083        saveModel=NULL;
1084      }
1085    }
1086    printf("Total time in crossover %g\n", CoinCpuTime()-startTime);
1087    delete [] saveUpper;
1088    delete [] saveLower;
1089  }
1090  return ;
1091}
1092#endif
1093/*****************************************************************************/
1094
1095// Default contructor
1096Idiot::Idiot()
1097{
1098  model_ = NULL;
1099  maxBigIts_ = 3;
1100  maxIts_ = 5;
1101  logLevel_ = 1; 
1102  logFreq_ = 100;
1103  maxIts2_ = 100;
1104  djTolerance_ = 1e-1;
1105  mu_ = 1e-4;
1106  drop_ = 5.0;
1107  exitDrop_=-1.0e20;
1108  muFactor_ = 0.3333;
1109  stopMu_ = 1e-12;
1110  smallInfeas_ = 1e-1;
1111  reasonableInfeas_ = 1e2;
1112  muAtExit_ =1.0e31;
1113  strategy_ =8;
1114  lambdaIterations_ =0;
1115  checkFrequency_ =100;
1116  whenUsed_ = NULL;
1117  majorIterations_ =30;
1118  exitFeasibility_ =-1.0;
1119  dropEnoughFeasibility_ =0.02;
1120  dropEnoughWeighted_ =0.01;
1121  // adjust
1122  double nrows=10000.0;
1123  int baseIts =(int) sqrt((double)nrows);
1124  baseIts =baseIts/10;
1125  baseIts *= 10;
1126  maxIts2_ =200+baseIts+5;
1127  reasonableInfeas_ =((double) nrows)*0.05;
1128  lightWeight_=0;
1129}
1130// Constructor from model
1131Idiot::Idiot(OsiSolverInterface &model)
1132{
1133  model_ = & model;
1134  maxBigIts_ = 3;
1135  maxIts_ = 5;
1136  logLevel_ = 1; 
1137  logFreq_ = 100;
1138  maxIts2_ = 100;
1139  djTolerance_ = 1e-1;
1140  mu_ = 1e-4;
1141  drop_ = 5.0;
1142  exitDrop_=-1.0e20;
1143  muFactor_ = 0.3333;
1144  stopMu_ = 1e-12;
1145  smallInfeas_ = 1e-1;
1146  reasonableInfeas_ = 1e2;
1147  muAtExit_ =1.0e31;
1148  strategy_ =8;
1149  lambdaIterations_ =0;
1150  checkFrequency_ =100;
1151  whenUsed_ = NULL;
1152  majorIterations_ =30;
1153  exitFeasibility_ =-1.0;
1154  dropEnoughFeasibility_ =0.02;
1155  dropEnoughWeighted_ =0.01;
1156  // adjust
1157  double nrows;
1158  if (model_)
1159    nrows=model_->getNumRows();
1160  else
1161    nrows=10000.0;
1162  int baseIts =(int) sqrt((double)nrows);
1163  baseIts =baseIts/10;
1164  baseIts *= 10;
1165  maxIts2_ =200+baseIts+5;
1166  reasonableInfeas_ =((double) nrows)*0.05;
1167  lightWeight_=0;
1168}
1169// Copy constructor.
1170Idiot::Idiot(const Idiot &rhs)
1171{
1172  model_ = rhs.model_;
1173  if (model_&&rhs.whenUsed_) {
1174    int numberColumns = model_->getNumCols();
1175    whenUsed_ = new int [numberColumns];
1176    memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1177  } else {
1178    whenUsed_=NULL;
1179  }
1180  djTolerance_ = rhs.djTolerance_;
1181  mu_ = rhs.mu_;
1182  drop_ = rhs.drop_;
1183  muFactor_ = rhs.muFactor_;
1184  stopMu_ = rhs.stopMu_;
1185  smallInfeas_ = rhs.smallInfeas_;
1186  reasonableInfeas_ = rhs.reasonableInfeas_;
1187  exitDrop_ = rhs.exitDrop_;
1188  muAtExit_ = rhs.muAtExit_;
1189  exitFeasibility_ = rhs.exitFeasibility_;
1190  dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1191  dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1192  maxBigIts_ = rhs.maxBigIts_;
1193  maxIts_ = rhs.maxIts_;
1194  majorIterations_ = rhs.majorIterations_;
1195  logLevel_ = rhs.logLevel_;
1196  logFreq_ = rhs.logFreq_;
1197  checkFrequency_ = rhs.checkFrequency_;
1198  lambdaIterations_ = rhs.lambdaIterations_;
1199  maxIts2_ = rhs.maxIts2_;
1200  strategy_ = rhs.strategy_;
1201  lightWeight_=rhs.lightWeight_;
1202}
1203// Assignment operator. This copies the data
1204Idiot & 
1205Idiot::operator=(const Idiot & rhs)
1206{
1207  if (this != &rhs) {
1208    delete [] whenUsed_;
1209    model_ = rhs.model_;
1210    if (model_&&rhs.whenUsed_) {
1211      int numberColumns = model_->getNumCols();
1212      whenUsed_ = new int [numberColumns];
1213      memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1214    } else {
1215      whenUsed_=NULL;
1216    }
1217    djTolerance_ = rhs.djTolerance_;
1218    mu_ = rhs.mu_;
1219    drop_ = rhs.drop_;
1220    muFactor_ = rhs.muFactor_;
1221    stopMu_ = rhs.stopMu_;
1222    smallInfeas_ = rhs.smallInfeas_;
1223    reasonableInfeas_ = rhs.reasonableInfeas_;
1224    exitDrop_ = rhs.exitDrop_;
1225    muAtExit_ = rhs.muAtExit_;
1226    exitFeasibility_ = rhs.exitFeasibility_;
1227    dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1228    dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1229    maxBigIts_ = rhs.maxBigIts_;
1230    maxIts_ = rhs.maxIts_;
1231    majorIterations_ = rhs.majorIterations_;
1232    logLevel_ = rhs.logLevel_;
1233    logFreq_ = rhs.logFreq_;
1234    checkFrequency_ = rhs.checkFrequency_;
1235    lambdaIterations_ = rhs.lambdaIterations_;
1236    maxIts2_ = rhs.maxIts2_;
1237    strategy_ = rhs.strategy_;
1238    lightWeight_=rhs.lightWeight_;
1239  }
1240  return *this;
1241}
1242Idiot::~Idiot()
1243{
1244  delete [] whenUsed_;
1245}
Note: See TracBrowser for help on using the repository browser.