source: trunk/Idiot.cpp @ 284

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

Fixes

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