source: branches/pre/Idiot.cpp @ 221

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

Still trying to go faster

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