source: branches/pre/Idiot.cpp @ 212

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

lots of stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.2 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  while ((mu>stopMu&&lastResult.infeas>smallInfeas)||
358         (lastResult.infeas<=smallInfeas&&
359         dropping(lastResult,exitDrop,smallInfeas,&badIts))||
360         checkIteration<2||lambdaIteration<lambdaIterations_) {
361    double keepinfeas;
362    if (lastResult.infeas<=exitFeasibility_)
363      break; 
364    iteration++;
365    checkIteration++;
366    if (lastResult.infeas<=smallInfeas&&lastResult.objval<bestFeasible) {
367      bestFeasible=lastResult.objval;
368    }
369    if ((saveStrategy&4096)) strategy |=256;
370    if ((saveStrategy&4)!=0&&iteration>2) {
371      /* go to relative tolerances */
372      double weighted=10.0+fabs(lastWeighted);
373      djExit=djTolerance_*weighted;
374      djFlag=2.0*djExit;
375      drop=drop_*weighted;
376      djTol=0.01*djExit;
377    }
378    result=IdiSolve(nrows,ncols,rowsol ,colsol,pi,dj,
379                     cost,rowlower,rowupper,
380                     lower,upper,elemXX,row,columnStart,columnLength,lambda,
381                     maxIts,mu,drop,
382                     maxmin,offset,strategy,djTol,djExit,djFlag);
383    n=0;
384    for (i=ordStart;i<ordEnd;i++) {
385      if (colsol[i]>lower[i]+fixTolerance) {
386        if (colsol[i]<upper[i]-fixTolerance) {
387          n++;
388        } else {
389          colsol[i]=upper[i];
390        }
391        whenUsed[i]=iteration;
392      } else {
393        colsol[i]=lower[i];
394      }
395    }
396    if ((logLevel_&1)!=0) {
397#ifdef CLP_IDIOT
398      if (!handler) {
399#endif
400        printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
401               iteration,result.infeas,result.objval,mu,result.iteration,n);
402#ifdef CLP_IDIOT
403      } else {
404        handler->message(CLP_IDIOT_ITERATION,*messages)
405          <<iteration<<result.infeas<<result.objval<<mu<<result.iteration<<n
406          <<CoinMessageEol;
407      }
408#endif
409    }
410    if (iteration>50&&n==numberAway&&result.infeas<1.0e-4)
411      break; // not much happening
412    numberAway=n;
413    keepinfeas = result.infeas;
414    lastWeighted=result.weighted;
415    iterationTotal += result.iteration;
416    if (iteration==1) {
417      if ((strategy_&1024)!=0&&mu<1.0e-10) 
418        result.infeas=firstInfeas*0.8;
419      if (result.infeas>firstInfeas*0.9
420          &&result.infeas>reasonableInfeas) {
421        iteration--;
422        memcpy(colsol,saveSol,ncols*sizeof(double));
423        mu*=1.0e-1;
424        bestFeasible=1.0e31;
425        bestWeighted=1.0e60;
426        if (mu<1.0e-30) break;
427      } else {
428        delete [] saveSol;
429        saveSol=0;
430        maxIts=maxIts2;
431        checkIteration=0;
432        if ((strategy_&1024)!=0) mu *= 1.0e-1;
433      }
434    }
435    if (iteration) {
436      /* this code is in to force it to terminate sometime */
437      double changeMu=factor;
438      if ((saveStrategy&64)!=0) {
439        keepinfeas=0.0; /* switch off ranga's increase */
440        fakeSmall=smallInfeas;
441      } else {
442        fakeSmall=-1.0;
443      }
444      saveLambdaScale=0.0;
445      if (result.infeas>reasonableInfeas||
446          (nTry+1==maxBigIts&&result.infeas>fakeSmall)) {
447        if (result.infeas>lastResult.infeas*(1.0-dropEnoughFeasibility_)||
448            nTry+1==maxBigIts||
449            (result.infeas>lastResult.infeas*0.9
450             &&result.weighted>lastResult.weighted
451             -dropEnoughWeighted_*fabs(lastResult.weighted))) {
452          mu*=changeMu;
453          if ((saveStrategy&32)!=0&&result.infeas<reasonableInfeas) {
454            reasonableInfeas=max(smallInfeas,reasonableInfeas*sqrt(changeMu));
455            printf("reasonable infeas now %g\n",reasonableInfeas);
456          }
457          nTry=0;
458          bestFeasible=1.0e31;
459          bestWeighted=1.0e60;
460          checkIteration=0;
461          lambdaIteration=0;
462#define LAMBDA
463#ifdef LAMBDA
464          if ((saveStrategy&2048)==0) {
465            memset(lambda,0,nrows*sizeof(double));
466          }
467#else
468          memset(lambda,0,nrows*sizeof(double));
469#endif
470        } else {
471          nTry++;
472        }
473      } else if (lambdaIterations_>=0) {
474        /* update lambda  */
475        double scale=1.0/mu;
476        int i,nnz=0;
477        saveLambdaScale=scale;
478         lambdaIteration++;
479         if ((saveStrategy&4)==0) drop = drop_/50.0;
480         if (lambdaIteration>4 && 
481            (((lambdaIteration%10)==0 && smallInfeas<keepinfeas) ||
482             (lambdaIteration%5)==0 && 1.5*smallInfeas<keepinfeas)) {
483           //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
484           smallInfeas *= 1.5;
485         }
486         if ((saveStrategy&2048)==0) {
487           for (i=0;i<nrows;i++) {
488             if (lambda[i]) nnz++;
489             lambda[i]+= scale*rowsol[i];
490           }
491         } else {
492           nnz=1;
493#ifdef LAMBDA
494           for (i=0;i<nrows;i++) {
495             lambda[i]+= scale*rowsol[i];
496           }
497#else
498           for (i=0;i<nrows;i++) {
499             lambda[i] = scale*rowsol[i];
500           }
501           for (i=0;i<ncols;i++) {
502             CoinBigIndex j;
503             double value=cost[i]*maxmin;
504             for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
505               int irow=row[j];
506               value+=element[j]*lambda[irow];
507             }
508             cost[i]=value*maxmin;
509           }
510           for (i=0;i<nrows;i++) {
511             offset+=lambda[i]*rowupper[i];
512             lambda[i]=0.0;
513           }
514#ifdef DEBUG
515           printf("offset %g\n",offset);
516#endif
517           model_->setDblParam(OsiObjOffset,offset);
518#endif
519         }
520        nTry++;
521        if (!nnz) {
522          bestFeasible=1.0e32;
523          bestWeighted=1.0e60;
524          checkIteration=0;
525          result.weighted=1.0e31;
526        }
527#ifdef DEBUG
528        for (i=0;i<ncols;i++) {
529          int j;
530          trueCost+=cost[i]*colsol[i];
531        }
532        printf("True objective %g\n",trueCost-offset);
533#endif
534      } else {
535        nTry++;
536      }
537      lastResult=result;
538      if (result.infeas<reasonableInfeas&&!belowReasonable) {
539        belowReasonable=1;
540        bestFeasible=1.0e32;
541        bestWeighted=1.0e60;
542        checkIteration=0;
543        result.weighted=1.0e31;
544      }
545    }
546    if (iteration>=majorIterations_) {
547      if (logLevel>2) 
548        printf("Exiting due to number of major iterations\n");
549      break;
550    }
551  }
552#define TRYTHIS
553#ifdef TRYTHIS
554  if ((saveStrategy&2048)!=0) {
555    double offset;
556    model_->getDblParam(OsiObjOffset,offset);
557    for (i=0;i<ncols;i++) {
558      CoinBigIndex j;
559      double djval=cost[i]*maxmin;
560      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
561        int irow=row[j];
562        djval -= element[j]*lambda[irow];
563      }
564      cost[i]=djval;
565    }
566    for (i=0;i<nrows;i++) {
567      offset+=lambda[i]*rowupper[i];
568    }
569    model_->setDblParam(OsiObjOffset,offset);
570  }
571#endif
572  if (saveLambdaScale) {
573    /* back off last update */
574    for (i=0;i<nrows;i++) {
575      lambda[i]-= saveLambdaScale*rowsol[i];
576    }
577  }
578  muAtExit_=mu;
579  n=0;
580  for (i=ordStart;i<ordEnd;i++) {
581    if (colsol[i]>lower[i]+fixTolerance) {
582      n++;
583      whenUsed[i]=iteration;
584    } else {
585      colsol[i]=lower[i];
586    }
587  }
588  if ((logLevel&1)==0) {
589    printf(
590            "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
591            iteration,mu,lastResult.infeas,lastResult.objval,n);
592  }
593  {
594    double large=0.0;
595    int i;
596    memset(rowsol,0,nrows*sizeof(double));
597    for (i=0;i<ncols;i++) {
598      CoinBigIndex j;
599      double value=colsol[i];
600      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
601        int irow=row[j];
602        rowsol[irow] += element[j]*value;
603      }
604    }
605    for (i=0;i<nrows;i++) {
606      if (rowsol[i] > rowupper[i]) {
607        double diff=rowsol[i]-rowupper[i];
608        if (diff>large) 
609          large=diff;
610      } else if (rowsol[i] < rowlower[i]) {
611        double diff=rowlower[i]-rowsol[i];
612        if (diff>large) 
613          large=diff;
614      } 
615    }
616    if (logLevel>2)
617      printf("largest infeasibility is %g\n",large);
618  }
619  /* subtract out lambda */
620  for (i=0;i<nrows;i++) {
621    pi[i]-=lambda[i];
622  }
623  for (i=0;i<ncols;i++) {
624    CoinBigIndex j;
625    double djval=cost[i]*maxmin;
626    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
627      int irow=row[j];
628      djval -= element[j]*pi[irow];
629    }
630    dj[i]=djval;
631  }
632  if ((strategy_&1024)!=0) {
633    double ratio = ((double) ncols)/((double) nrows);
634    printf("col/row ratio %g infeas ratio %g\n",ratio,lastResult.infeas/firstInfeas);
635    if (lastResult.infeas>0.01*firstInfeas*ratio) {
636      strategy_ &= (~1024);
637      printf(" - layer off\n");
638    } else {
639      printf(" - layer on\n");
640    }
641  }
642  delete [] saveSol;
643  delete [] lambda;
644  // save solution
645  // duals not much use - but save anyway
646#ifdef CLP_IDIOT
647  memcpy(model_->primalRowSolution(),rowsol,nrows*sizeof(double));
648  memcpy(model_->primalColumnSolution(),colsol,ncols*sizeof(double));
649  memcpy(model_->dualRowSolution(),pi,nrows*sizeof(double));
650  memcpy(model_->dualColumnSolution(),dj,ncols*sizeof(double));
651#else
652  model_->setColSolution(colsol);
653  model_->setRowPrice(pi);
654  delete [] cost;
655#endif
656  delete [] rowsol;
657  delete [] colsol;
658  delete [] pi;
659  delete [] dj;
660  return ;
661}
662#ifdef CLP_IDIOT
663void
664Idiot::crossOver(int mode)
665{
666  double fixTolerance=IDIOT_FIX_TOLERANCE;
667  double startTime = CoinCpuTime();
668  ClpSimplex * saveModel=NULL;
669  ClpMatrixBase * matrix = model_->clpMatrix();
670  const int * row = matrix->getIndices();
671  const CoinBigIndex * columnStart = matrix->getVectorStarts();
672  const int * columnLength = matrix->getVectorLengths(); 
673  const double * element = matrix->getElements();
674  const double * rowupper = model_->getRowUpper();
675  int nrows=model_->getNumRows();
676  int ncols=model_->getNumCols();
677  double * rowsol, * colsol;
678  // different for Osi
679  double * lower = model_->columnLower();
680  double * upper = model_->columnUpper();
681  const double * rowlower= model_->getRowLower();
682  int * whenUsed=whenUsed_;
683  rowsol= model_->primalRowSolution();
684  colsol= model_->primalColumnSolution();;
685  double * cost=model_->objective();
686
687  int slackEnd,ordStart,ordEnd;
688  int slackStart = countCostedSlacks(model_);
689
690  int addAll = mode&7;
691  int presolve=0;
692
693  double djTolerance = djTolerance_;
694  if (djTolerance>0.0&&djTolerance<1.0)
695    djTolerance=1.0;
696  int iteration;
697  int i, n=0;
698  double ratio=1.0;
699  double objValue=0.0;
700  if ((strategy_&128)!=0) {
701    fixTolerance=SMALL_IDIOT_FIX_TOLERANCE;
702  }
703  if ((mode&16)!=0&&addAll<2) presolve=1;
704  double * saveUpper = NULL;
705  double * saveLower = NULL;
706  if (addAll<3) {
707    saveUpper = new double [ncols];
708    saveLower = new double [ncols];
709    memcpy(saveUpper,upper,ncols*sizeof(double));
710    memcpy(saveLower,lower,ncols*sizeof(double));
711  }
712  if (slackStart>=0) {
713    slackEnd=slackStart+nrows;
714    if (slackStart) {
715      ordStart=0;
716      ordEnd=slackStart;
717    } else {
718      ordStart=nrows;
719      ordEnd=ncols;
720    }
721  } else {
722    slackEnd=slackStart;
723    ordStart=0;
724    ordEnd=ncols;
725  }
726  /* get correct rowsol (without known slacks) */
727  memset(rowsol,0,nrows*sizeof(double));
728  for (i=ordStart;i<ordEnd;i++) {
729    CoinBigIndex j;
730    double value=colsol[i];
731    if (value<lower[i]+fixTolerance) {
732      value=lower[i];
733      colsol[i]=value;
734    }
735    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
736      int irow=row[j];
737      rowsol[irow]+=value*element[j];
738    }
739  }
740  if (slackStart>=0) {
741    for (i=0;i<nrows;i++) {
742      if (ratio*rowsol[i]>rowlower[i]&&rowsol[i]>1.0e-8) {
743        ratio=rowlower[i]/rowsol[i];
744      }
745    }
746    for (i=0;i<nrows;i++) {
747      rowsol[i]*=ratio;
748    }
749    for (i=ordStart;i<ordEnd;i++) {
750      double value=colsol[i]*ratio;
751      colsol[i]=value;
752      objValue+=value*cost[i];
753    }
754    for (i=0;i<nrows;i++) {
755      double value=rowlower[i]-rowsol[i];
756      colsol[i+slackStart]=value;
757      objValue+=value*cost[i+slackStart];
758    }
759    printf("New objective after scaling %g\n",objValue);
760  }
761  model_->createStatus();
762  /* addAll
763     0 - chosen,all used, all
764     1 - chosen, all
765     2 - all
766     3 - do not do anything  - maybe basis
767  */
768  for (i=ordStart;i<ordEnd;i++) {
769    if (addAll<2) {
770      if (colsol[i]<lower[i]+fixTolerance) {
771        upper[i]=lower[i];
772        colsol[i]=lower[i];
773      } else if (colsol[i]>upper[i]-fixTolerance) {
774        lower[i]=upper[i];
775        colsol[i]=upper[i];
776      }
777    }
778    model_->setColumnStatus(i,ClpSimplex::superBasic);
779  }
780  if (slackStart>=0) {
781    for (i=0;i<nrows;i++) {
782      model_->setRowStatus(i,ClpSimplex::superBasic);
783    }
784    for (i=slackStart;i<slackEnd;i++) {
785      model_->setColumnStatus(i,ClpSimplex::basic);
786    }
787  } else {
788    /* still try and put singletons rather than artificials in basis */
789    int ninbas=0;
790    for (i=0;i<nrows;i++) {
791      model_->setRowStatus(i,ClpSimplex::basic);
792    }
793    for (i=0;i<ncols;i++) {
794      if (columnLength[i]==1&&upper[i]>lower[i]+1.0e-5) {
795        CoinBigIndex j =columnStart[i];
796        double value=element[j];
797        int irow=row[j];
798        double rlo=rowlower[irow];
799        double rup=rowupper[irow];
800        double clo=lower[i];
801        double cup=upper[i];
802        double csol=colsol[i];
803        /* adjust towards feasibility */
804        double move=0.0;
805        if (rowsol[irow]>rup) {
806          move=(rup-rowsol[irow])/value;
807          if (value>0.0) {
808            /* reduce */
809            if (csol+move<clo) move=min(0.0,clo-csol);
810          } else {
811            /* increase */
812            if (csol+move>cup) move=max(0.0,cup-csol);
813          }
814        } else if (rowsol[irow]<rlo) {
815          move=(rlo-rowsol[irow])/value;
816          if (value>0.0) {
817            /* increase */
818            if (csol+move>cup) move=max(0.0,cup-csol);
819          } else {
820            /* reduce */
821            if (csol+move<clo) move=min(0.0,clo-csol);
822          }
823        } else {
824          /* move to improve objective */
825          if (cost[i]>=0.0) {
826            if (value>0.0) {
827              move=(rlo-rowsol[irow])/value;
828              /* reduce */
829              if (csol+move<clo) move=min(0.0,clo-csol);
830            } else {
831              move=(rup-rowsol[irow])/value;
832              /* increase */
833              if (csol+move>cup) move=max(0.0,cup-csol);
834            }
835          } else {
836            if (value>0.0) {
837              move=(rup-rowsol[irow])/value;
838              /* increase */
839              if (csol+move>cup) move=max(0.0,cup-csol);
840            } else {
841              move=(rlo-rowsol[irow])/value;
842              /* reduce */
843              if (csol+move<clo) move=min(0.0,clo-csol);
844            }
845          }
846        }
847        rowsol[irow] +=move*value;
848        colsol[i]+=move;
849        /* put in basis if row was artificial */
850        if (rup-rlo<1.0e-7&&model_->getRowStatus(irow)==ClpSimplex::basic) {
851          model_->setRowStatus(irow,ClpSimplex::superBasic);
852          model_->setColumnStatus(i,ClpSimplex::basic);
853          ninbas++;
854        }
855      }
856    }
857    /*printf("%d in basis\n",ninbas);*/
858  }
859  if (addAll<3) {
860    ClpPresolve pinfo;
861    if (presolve) {
862      saveModel = model_;
863      model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
864    }
865    model_->primal(1);
866    if (presolve) {
867      pinfo.postsolve(true);
868      delete model_;
869      model_ = saveModel;
870      saveModel=NULL;
871    }
872    if (addAll<2) {
873      n=0;
874      if (!addAll ) {
875        /* could do scans to get a good number */
876        iteration=1;
877        for (i=ordStart;i<ordEnd;i++) {
878          if (whenUsed[i]>=iteration) {
879            if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
880              n++;
881              upper[i]=saveUpper[i];
882              lower[i]=saveLower[i];
883            }
884          }
885        }
886      } else {
887        for (i=ordStart;i<ordEnd;i++) {
888          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
889            n++;
890            upper[i]=saveUpper[i];
891            lower[i]=saveLower[i];
892          }
893        }
894      }
895      printf("Time so far %g, %d now added from previous iterations\n",
896             CoinCpuTime()-startTime,n);
897      if (presolve) {
898        saveModel = model_;
899        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
900      } else {
901        presolve=0;
902      }
903      model_->primal(1);
904      if (presolve) {
905        pinfo.postsolve(true);
906        delete model_;
907        model_ = saveModel;
908        saveModel=NULL;
909      }
910      if (!addAll) {
911        n=0;
912        for (i=ordStart;i<ordEnd;i++) {
913          if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
914            n++;
915            upper[i]=saveUpper[i];
916            lower[i]=saveLower[i];
917          }
918        }
919        printf("Time so far %g, %d now added from previous iterations\n",
920               CoinCpuTime()-startTime,n);
921      }
922      if (presolve) {
923        saveModel = model_;
924        model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
925      } else {
926        presolve=0;
927      }
928      model_->primal(1);
929      if (presolve) {
930        pinfo.postsolve(true);
931        delete model_;
932        model_ = saveModel;
933        saveModel=NULL;
934      }
935    }
936    printf("Total time in crossover %g\n", CoinCpuTime()-startTime);
937    delete [] saveUpper;
938    delete [] saveLower;
939  }
940  return ;
941}
942#endif
943/*****************************************************************************/
944
945// Default contructor
946Idiot::Idiot()
947{
948  model_ = NULL;
949  maxBigIts_ = 3;
950  maxIts_ = 5;
951  logLevel_ = 1; 
952  logFreq_ = 100;
953  maxIts2_ = 100;
954  djTolerance_ = 1e-1;
955  mu_ = 1e-4;
956  drop_ = 5.0;
957  exitDrop_=-1.0e20;
958  muFactor_ = 0.3333;
959  stopMu_ = 1e-12;
960  smallInfeas_ = 1e-1;
961  reasonableInfeas_ = 1e2;
962  muAtExit_ =1.0e31;
963  strategy_ =8;
964  lambdaIterations_ =0;
965  checkFrequency_ =100;
966  whenUsed_ = NULL;
967  majorIterations_ =30;
968  exitFeasibility_ =-1.0;
969  dropEnoughFeasibility_ =0.02;
970  dropEnoughWeighted_ =0.01;
971  // adjust
972  double nrows=10000.0;
973  int baseIts =(int) sqrt((double)nrows);
974  baseIts =baseIts/10;
975  baseIts *= 10;
976  maxIts2_ =200+baseIts+5;
977  reasonableInfeas_ =((double) nrows)*0.05;
978}
979// Constructor from model
980Idiot::Idiot(OsiSolverInterface &model)
981{
982  model_ = & model;
983  maxBigIts_ = 3;
984  maxIts_ = 5;
985  logLevel_ = 1; 
986  logFreq_ = 100;
987  maxIts2_ = 100;
988  djTolerance_ = 1e-1;
989  mu_ = 1e-4;
990  drop_ = 5.0;
991  exitDrop_=-1.0e20;
992  muFactor_ = 0.3333;
993  stopMu_ = 1e-12;
994  smallInfeas_ = 1e-1;
995  reasonableInfeas_ = 1e2;
996  muAtExit_ =1.0e31;
997  strategy_ =8;
998  lambdaIterations_ =0;
999  checkFrequency_ =100;
1000  whenUsed_ = NULL;
1001  majorIterations_ =30;
1002  exitFeasibility_ =-1.0;
1003  dropEnoughFeasibility_ =0.02;
1004  dropEnoughWeighted_ =0.01;
1005  // adjust
1006  double nrows;
1007  if (model_)
1008    nrows=model_->getNumRows();
1009  else
1010    nrows=10000.0;
1011  int baseIts =(int) sqrt((double)nrows);
1012  baseIts =baseIts/10;
1013  baseIts *= 10;
1014  maxIts2_ =200+baseIts+5;
1015  reasonableInfeas_ =((double) nrows)*0.05;
1016}
1017// Copy constructor.
1018Idiot::Idiot(const Idiot &rhs)
1019{
1020  model_ = rhs.model_;
1021  if (model_&&rhs.whenUsed_) {
1022    int numberColumns = model_->getNumCols();
1023    whenUsed_ = new int [numberColumns];
1024    memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1025  } else {
1026    whenUsed_=NULL;
1027  }
1028  djTolerance_ = rhs.djTolerance_;
1029  mu_ = rhs.mu_;
1030  drop_ = rhs.drop_;
1031  muFactor_ = rhs.muFactor_;
1032  stopMu_ = rhs.stopMu_;
1033  smallInfeas_ = rhs.smallInfeas_;
1034  reasonableInfeas_ = rhs.reasonableInfeas_;
1035  exitDrop_ = rhs.exitDrop_;
1036  muAtExit_ = rhs.muAtExit_;
1037  exitFeasibility_ = rhs.exitFeasibility_;
1038  dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1039  dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1040  maxBigIts_ = rhs.maxBigIts_;
1041  maxIts_ = rhs.maxIts_;
1042  majorIterations_ = rhs.majorIterations_;
1043  logLevel_ = rhs.logLevel_;
1044  logFreq_ = rhs.logFreq_;
1045  checkFrequency_ = rhs.checkFrequency_;
1046  lambdaIterations_ = rhs.lambdaIterations_;
1047  maxIts2_ = rhs.maxIts2_;
1048  strategy_ = rhs.strategy_;
1049}
1050// Assignment operator. This copies the data
1051Idiot & 
1052Idiot::operator=(const Idiot & rhs)
1053{
1054  if (this != &rhs) {
1055    delete [] whenUsed_;
1056    model_ = rhs.model_;
1057    if (model_&&rhs.whenUsed_) {
1058      int numberColumns = model_->getNumCols();
1059      whenUsed_ = new int [numberColumns];
1060      memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
1061    } else {
1062      whenUsed_=NULL;
1063    }
1064    djTolerance_ = rhs.djTolerance_;
1065    mu_ = rhs.mu_;
1066    drop_ = rhs.drop_;
1067    muFactor_ = rhs.muFactor_;
1068    stopMu_ = rhs.stopMu_;
1069    smallInfeas_ = rhs.smallInfeas_;
1070    reasonableInfeas_ = rhs.reasonableInfeas_;
1071    exitDrop_ = rhs.exitDrop_;
1072    muAtExit_ = rhs.muAtExit_;
1073    exitFeasibility_ = rhs.exitFeasibility_;
1074    dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1075    dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1076    maxBigIts_ = rhs.maxBigIts_;
1077    maxIts_ = rhs.maxIts_;
1078    majorIterations_ = rhs.majorIterations_;
1079    logLevel_ = rhs.logLevel_;
1080    logFreq_ = rhs.logFreq_;
1081    checkFrequency_ = rhs.checkFrequency_;
1082    lambdaIterations_ = rhs.lambdaIterations_;
1083    maxIts2_ = rhs.maxIts2_;
1084    strategy_ = rhs.strategy_;
1085  }
1086  return *this;
1087}
1088Idiot::~Idiot()
1089{
1090  delete [] whenUsed_;
1091}
Note: See TracBrowser for help on using the repository browser.