source: branches/pre/Idiot.cpp @ 210

Last change on this file since 210 was 210, checked in by forrest, 17 years ago

Trying

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