source: trunk/Idiot.cpp @ 76

Last change on this file since 76 was 76, checked in by jpfasano, 17 years ago

Fix to compile on windows

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