source: trunk/Idiot.cpp @ 403

Last change on this file since 403 was 403, checked in by forrest, 16 years ago

to stop Idiot seg fault

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