source: branches/pre/Idiot.cpp @ 1926

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

Stuff

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