source: trunk/IdiSolve.cpp @ 225

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

This should break everything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.4 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 "CoinHelperFunctions.hpp"
10#include "Idiot.hpp"
11#define FIT
12#ifdef FIT
13#define HISTORY 8
14#else
15#define HISTORY 7
16#endif
17#define NSOLVE HISTORY-1
18static void solveSmall(int nsolve,double **aIn,double **a, double * b) {
19  int i,j;
20  /* copy */
21  for (i=0;i<nsolve;i++) {
22    for (j=0;j<nsolve;j++) {
23      a[i][j]=aIn[i][j];
24    }
25  }
26  for (i=0;i<nsolve;i++) {
27    /* update using all previous */
28    double diagonal;
29    int j;
30    for (j=i;j<nsolve;j++) {
31      int k;
32      double value=a[i][j];
33      for (k=0;k<i;k++) {
34        value-=a[k][i]*a[k][j];
35      }
36      a[i][j]=value;
37    }
38    diagonal=a[i][i];
39    if (diagonal<1.0e-20) {
40      diagonal=0.0;
41    } else {
42      diagonal=1.0/sqrt(diagonal);
43    }
44    a[i][i]=diagonal;
45    for (j=i+1;j<nsolve;j++) {
46      a[i][j]*=diagonal;
47    }
48  }
49  /* update */
50  for (i=0;i<nsolve;i++) {
51    int j;
52    double value=b[i];
53    for (j=0;j<i;j++) {
54      value-=b[j]*a[j][i];
55    }
56    value*=a[i][i];
57    b[i]=value;
58  }
59  for (i=nsolve-1;i>=0;i--) {
60    int j;
61    double value=b[i];
62    for (j=i+1;j<nsolve;j++) {
63      value-=b[j]*a[i][j];
64    }
65    value*=a[i][i];
66    b[i]=value;
67  }
68}
69IdiotResult
70Idiot::objval(int nrows, int ncols, double * rowsol , double * colsol,
71              double * pi, double * djs, const double * cost , 
72                        const double * rowlower,
73              const double * rowupper, const double * lower,
74              const double * upper, const double * elemnt, 
75              const int * row, const CoinBigIndex * columnStart,
76                          const int * length, int extraBlock, int * rowExtra,
77                 double * solExtra, double * elemExtra, double * upperExtra,
78                 double * costExtra,double weight)
79{
80  IdiotResult result;
81  double objvalue=0.0;
82  double sum1=0.0,sum2=0.0;
83  int i;
84  for (i=0;i<nrows;i++) {
85    rowsol[i]=-rowupper[i];
86  }
87  for (i=0;i<ncols;i++) {
88    CoinBigIndex j;
89    double value=colsol[i];
90    if (value) {
91      objvalue+=value*cost[i];
92      if (elemnt) {
93        for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
94          int irow=row[j];
95          rowsol[irow]+=elemnt[j]*value;
96        }
97      } else {
98        for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
99          int irow=row[j];
100          rowsol[irow]+=value;
101        }
102      }
103    }
104  }
105  /* adjust to make as feasible as possible */
106  /* no */
107  if (extraBlock) {
108    for (i=0;i<extraBlock;i++) {
109      double element=elemExtra[i];
110      int irow=rowExtra[i];
111      objvalue+=solExtra[i]*costExtra[i];
112      rowsol[irow] += solExtra[i]*element;
113    }
114  }
115  for (i=0;i<nrows;i++) {
116    double value=rowsol[i];
117    sum1+=fabs(value);
118    sum2+=value*value;
119    pi[i]=-2.0*weight*value;
120  }
121  result.infeas=sum1;
122  result.objval=objvalue;
123  result.weighted=objvalue+weight*sum2;
124  result.sumSquared=sum2;
125  return result;
126}
127IdiotResult
128Idiot::IdiSolve(
129                    int nrows, int ncols, double * rowsol , double * colsol,
130              double * pi, double * djs, const double * origcost , const double * rowlower,
131              double * rowupper, const double * lower,
132              const double * upper, const double * elemnt, 
133                    const int * row, const CoinBigIndex * columnStart,
134                    const int * length, double * lambda,
135                    int maxIts,double mu,double drop,
136                    double maxmin, double offset,
137                    int strategy,double djTol,double djExit,double djFlag)
138{
139  IdiotResult result;
140  int  i,j,k,iter;
141  double value=0.0,objvalue=0.0,weightedObj=0.0;
142  double tolerance=1.0e-8;
143  double * history[HISTORY+1];
144  int ncolx;
145  int nChange;
146  int extraBlock=0;
147  int * rowExtra=NULL;
148  double * solExtra=NULL;
149  double * elemExtra=NULL;
150  double * upperExtra=NULL;
151  double * costExtra=NULL;
152  double * useCostExtra=NULL;
153  double * saveExtra=NULL;
154  double * cost=NULL;
155  double saveValue=1.0e30;
156  double saveOffset=offset;
157  double useOffset=offset;
158  /*#define NULLVECTOR*/
159#ifndef NULLVECTOR
160  int nsolve=NSOLVE;
161#else
162  int nsolve=NSOLVE+1;/* allow for null vector */
163#endif
164  int nflagged;
165  double * thetaX;
166  double * djX;
167  double * bX;
168  double * vX;
169  double ** aX;
170  double **aworkX;
171  double ** allsum;
172  double * saveSol=0;
173  const double * useCost=cost;
174  double bestSol=1.0e60;
175  double weight=0.5/mu;
176  char * statusSave=new char[2*ncols];
177  char * statusWork=statusSave+ncols;
178#define DJTEST 5
179  double djSave[DJTEST];
180  double largestDj=0.0;
181  double smallestDj=1.0e60;
182  double maxDj=0.0;
183  int doFull=0;
184#define SAVEHISTORY 10
185#define EVERY (2*SAVEHISTORY)
186#define AFTER SAVEHISTORY*(HISTORY+1)
187#define DROP 5
188  double after=AFTER;
189  double obj[DROP];
190  double kbad=0,kgood=0;
191  if (strategy&128) after=999999; /* no acceleration at all */
192  for (i=0;i<DROP;i++) {
193    obj[i]=1.0e70;
194  }
195  allsum=new double * [nsolve];
196  aX=new double * [nsolve];
197  aworkX=new double * [nsolve];
198  thetaX=new double[nsolve];
199  vX=new double[nsolve];
200  bX=new double[nsolve];
201  djX=new double[nsolve];
202  allsum[0]=pi;
203  for (i=0;i<nsolve;i++) {
204    if (i) allsum[i]=new double[nrows];
205    aX[i]=new double[nsolve];
206    aworkX[i]=new double[nsolve];
207  }
208  /* check = rows */
209  for (i=0;i<nrows;i++) {
210    if (rowupper[i]-rowlower[i]>tolerance) {
211      extraBlock++;
212    }
213  }
214  cost=new double[ncols];
215  memset(rowsol,0,nrows*sizeof(double));
216  for (i=0;i<ncols;i++) {
217    CoinBigIndex j;
218    double value=origcost[i]*maxmin;
219    double value2=colsol[i];
220    if (elemnt) {
221      for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
222        int irow=row[j];
223        value+=elemnt[j]*lambda[irow];
224        rowsol[irow]+=elemnt[j]*value2;
225      }
226    } else {
227      for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
228        int irow=row[j];
229        value+=lambda[irow];
230        rowsol[irow]+=value2;
231      }
232    }
233    cost[i]=value;
234  }
235  if (extraBlock) {
236    rowExtra=new int[extraBlock];
237    solExtra=new double[extraBlock];
238    elemExtra=new double[extraBlock];
239    upperExtra=new double[extraBlock];
240    costExtra=new double[extraBlock];
241    saveExtra=new double[extraBlock];
242    extraBlock=0;
243    for (i=0;i<nrows;i++) {
244      if (rowupper[i]-rowlower[i]>tolerance) {
245        double smaller,difference;
246        double value;
247        saveExtra[extraBlock]=rowupper[i];
248        if (fabs(rowupper[i])>fabs(rowlower[i])) {
249          smaller = rowlower[i];
250          value=-1.0;
251        } else {
252          smaller = rowupper[i];
253          value=1.0;
254        }
255        if (fabs(smaller)>1.0e10) {
256          printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
257                 i,smaller);
258          abort();
259        }
260        difference=rowupper[i]-rowlower[i];
261        difference = min(difference,1.0e31);
262        rowupper[i]=smaller;
263        elemExtra[extraBlock]=value;
264        solExtra[extraBlock]=(rowupper[i]-rowsol[i])/value;
265        if (solExtra[extraBlock]<0.0) solExtra[extraBlock]=0.0;
266        if (solExtra[extraBlock]>difference) solExtra[extraBlock]=difference;
267        costExtra[extraBlock]=lambda[i]*value;
268        upperExtra[extraBlock]=difference;
269        rowsol[i]+=value*solExtra[extraBlock];
270        rowExtra[extraBlock++]=i;
271      }
272    }
273  }
274  for (i=0;i<nrows;i++) {
275    offset+=lambda[i]*rowsol[i];
276  }
277  if ((strategy&256)!=0) {
278    /* save best solution */
279    saveSol=new double[ncols];
280    memcpy(saveSol,colsol,ncols*sizeof(double));
281    if (extraBlock) {
282      useCostExtra=new double[extraBlock];
283      memset(useCostExtra,0,extraBlock*sizeof(double));
284    }
285    useCost=origcost;
286    useOffset=saveOffset;
287  } else {
288    useCostExtra=costExtra;
289    useCost=cost;
290    useOffset=offset;
291  }
292  ncolx=ncols+extraBlock;
293  for (i=0;i<HISTORY+1;i++) {
294    history[i]=new double[ncolx];
295  }
296  for (i=0;i<DJTEST;i++) {
297    djSave[i]=1.0e30;
298  }
299  for (i=0;i<ncols;i++) {
300    if (upper[i]-lower[i]) {
301      statusSave[i]=0;
302    } else {
303      statusSave[i]=1;
304    }
305  }
306  // for two pass method
307  int start[2];
308  int stop[2];
309  int direction=-1;
310  start[0]=0;stop[0]=ncols;
311  start[1]=0;stop[1]=0;
312  iter=0;
313  for (;iter<maxIts;iter++) {
314    double sum1=0.0,sum2=0.0,djtot;
315    double lastObj=1.0e70;
316    int good=0,doScale=0;
317    if (strategy&16) {
318      int ii=iter/EVERY+1;
319      ii=ii*EVERY;
320      if (iter>ii-HISTORY*2&&(iter&1)==0) {
321        double * x=history[HISTORY-1];
322        for (i=HISTORY-1;i>0;i--) {
323          history[i]=history[i-1];
324        }
325        history[0]=x;
326        memcpy(history[0],colsol,ncols*sizeof(double));
327        memcpy(history[0]+ncols,solExtra,extraBlock*sizeof(double));
328      }
329    }
330    if ((iter % SAVEHISTORY)==0||doFull) {
331      if ((strategy&16)==0) {
332        double * x=history[HISTORY-1];
333        for (i=HISTORY-1;i>0;i--) {
334          history[i]=history[i-1];
335        }
336        history[0]=x;
337        memcpy(history[0],colsol,ncols*sizeof(double));
338        memcpy(history[0]+ncols,solExtra,extraBlock*sizeof(double));
339      }
340    }
341    /* start full try */
342    if ((iter % EVERY)==0||doFull) {
343      // for next pass
344      direction = - direction;
345      // randomize.
346      // The cast is to avoid gcc compiler warning
347      int kcol = (int)(ncols*CoinDrand48());
348      if (kcol==ncols)
349        kcol=ncols-1;
350      if (direction>0) {
351        start[0]=kcol;stop[0]=ncols;
352        start[1]=0;stop[1]=kcol;
353      } else {
354        start[0]=kcol;stop[0]=-1;
355        start[1]=ncols-1;stop[1]=kcol;
356      }
357      int itry=0;
358      /*if ((strategy&16)==0) {
359        double * x=history[HISTORY-1];
360        for (i=HISTORY-1;i>0;i--) {
361          history[i]=history[i-1];
362        }
363        history[0]=x;
364        memcpy(history[0],colsol,ncols*sizeof(double));
365        memcpy(history[0]+ncols,solExtra,extraBlock*sizeof(double));
366        }*/
367      while (!good) {
368        itry++;
369#define MAXTRY 5
370        if (iter>after&&doScale<2&&itry<MAXTRY) {
371          /* now full one */
372          for (i=0;i<nrows;i++) {
373            rowsol[i]=-rowupper[i];
374          }
375          sum2=0.0;
376          objvalue=0.0;
377          memset(pi,0,nrows*sizeof(double));
378          djtot=0.0;
379          {
380            double * theta=thetaX;
381            double * dj=djX;
382            double * b=bX;
383            double ** a=aX;
384            double ** awork=aworkX;
385            double * v=vX;
386            double c;
387#ifdef FIT
388            int ntot=0,nsign=0,ngood=0,mgood[4]={0,0,0,0};
389            double diff1,diff2,val0,val1,val2,newValue;
390            memcpy(history[HISTORY-1],colsol,ncols*sizeof(double));
391            memcpy(history[HISTORY-1]+ncols,solExtra,extraBlock*sizeof(double));
392#endif
393            dj[0]=0.0;
394            for (i=1;i<nsolve;i++) {
395              dj[i]=0.0;
396              memset(allsum[i],0,nrows*sizeof(double));
397            }
398            for (i=0;i<ncols;i++) {
399              double value2=colsol[i];
400              if (value2>lower[i]+tolerance) {
401                if(value2<(upper[i]-tolerance)) {
402                  int k;
403                  objvalue+=value2*cost[i];
404#ifdef FIT
405                  ntot++;
406                  val0=history[0][i];
407                  val1=history[1][i];
408                  val2=history[2][i];
409                  diff1=val0-val1;
410                  diff2=val1-val2;
411                  if (diff1*diff2>=0.0) {
412                    nsign++;
413                    if (fabs(diff1)<=fabs(diff2)) {
414          // cast is to avoid gcc compiler
415          // warning: initialization to 'int' from 'double'
416                      int ii=(int)fabs(4.0*diff1/diff2);
417                      if (ii==4) ii=3;
418                      mgood[ii]++;
419                      ngood++;
420                    }
421                    if (fabs(diff1)<0.75*fabs(diff2)) {
422                      newValue=val1+(diff1*diff2)/(diff2-diff1);
423                    } else {
424                      newValue=val1+4.0*diff1;
425                    }
426                  } else {
427                    newValue=0.333333333*(val0+val1+val2);
428                  }
429                  if (newValue>upper[i]-tolerance) {
430                    newValue=upper[i];
431                  } else if (newValue<lower[i]+tolerance) {
432                    newValue=lower[i];
433                  }
434                  history[HISTORY-1][i]=newValue;
435#endif
436                  for (k=0;k<HISTORY-1;k++) {
437                    value=history[k][i]-history[k+1][i];
438                    dj[k] += value*cost[i];
439                    v[k]=value;
440                  }
441                  if (elemnt) {
442                    for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
443                      int irow=row[j];
444                      for (k=0;k<HISTORY-1;k++) {
445                        allsum[k][irow] += elemnt[j]*v[k];
446                      }
447                      rowsol[irow] += elemnt[j]*value2;
448                    }
449                  } else {
450                    for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
451                      int irow=row[j];
452                      for (k=0;k<HISTORY-1;k++) {
453                        allsum[k][irow] += v[k];
454                      }
455                      rowsol[irow] += value2;
456                    }
457                  }
458                } else {
459                  /* at ub */
460                  colsol[i]=upper[i];
461                  value2=colsol[i];
462                  objvalue+=value2*cost[i];
463                  if (elemnt) {
464                    for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
465                      int irow=row[j];
466                      rowsol[irow] += elemnt[j]*value2;
467                    }
468                  } else {
469                    for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
470                      int irow=row[j];
471                      rowsol[irow] += value2;
472                    }
473                  }
474                }
475              } else {
476                /* at lb */
477                if (value2) {
478                  objvalue+=value2*cost[i];
479                  if (elemnt) {
480                    for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
481                      int irow=row[j];
482                      rowsol[irow] += elemnt[j]*value2;
483                    }
484                  } else {
485                    for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
486                      int irow=row[j];
487                      rowsol[irow] += value2;
488                    }
489                  }
490                }
491              }
492            }
493#ifdef FIT
494            /*printf("total %d, same sign %d, good %d %d %d %d %d\n",
495              ntot,nsign,ngood,mgood[0],mgood[1],mgood[2],mgood[3]);*/
496#endif
497            if (extraBlock) {
498              for (i=0;i<extraBlock;i++) {
499                double element=elemExtra[i];
500                int irow=rowExtra[i];
501                objvalue+=solExtra[i]*costExtra[i];
502                if (solExtra[i]>tolerance
503                    &&solExtra[i]<(upperExtra[i]-tolerance)) {
504                  double value2=solExtra[i];
505                  int k;
506                  for (k=0;k<HISTORY-1;k++) {
507                    value=history[k][i+ncols]-history[k+1][i+ncols];
508                    dj[k] += value*costExtra[i];
509                    allsum[k][irow] += element*value;
510                  }
511                  rowsol[irow] += element*value2;
512                } else {
513                  double value2=solExtra[i];
514                  double element=elemExtra[i];
515                  int irow=rowExtra[i];
516                  rowsol[irow] += element*value2;
517                }
518              }
519            }
520#ifdef NULLVECTOR
521            if ((strategy&64)) {
522              double djVal=dj[0];
523              for (i=0;i<ncols-nrows;i++) {
524                double value2=colsol[i];
525                if (value2>lower[i]+tolerance&&value2<upper[i]-tolerance) {
526                  value=history[0][i]-history[1][i];
527                } else {
528                  value=0.0;
529                }
530                history[HISTORY][i]=value;
531              }
532              for (;i<ncols;i++) {
533                double value2=colsol[i];
534                double delta;
535                int irow=i-(ncols-nrows);
536                double oldSum=allsum[0][irow];;
537                if (value2>lower[i]+tolerance&&value2<upper[i]-tolerance) {
538                  delta=history[0][i]-history[1][i];
539                } else {
540                  delta=0.0;
541                }
542                djVal -= delta*cost[i];
543                oldSum -= delta;
544                delta = - oldSum;
545                djVal += delta*cost[i];
546                history[HISTORY][i]=delta;
547              }
548              dj[HISTORY-1]=djVal;
549              djVal=0.0;
550              for (i=0;i<ncols;i++) {
551                double value2=colsol[i];
552                if (value2>lower[i]+tolerance&&value2<upper[i]-tolerance||
553                    i>=ncols-nrows) {
554                  int k;
555                  value=history[HISTORY][i];
556                  djVal += value*cost[i];
557                  for (j=columnStart[i];j<columnStart[i]+length[i];j++) {
558                    int irow=row[j];
559                    allsum[nsolve-1][irow] += value;
560                  }
561                }
562              }
563              printf("djs %g %g\n",dj[HISTORY-1],djVal);
564            }
565#endif
566            for (i=0;i<nsolve;i++) {
567              int j;
568              b[i]=0.0;
569              for (j=0;j<nsolve;j++) {
570                a[i][j]=0.0;
571              }
572            }
573            c=0.0;     
574            for (i=0;i<nrows;i++) {
575              double value=rowsol[i];
576              for (k=0;k<nsolve;k++) {
577                v[k]=allsum[k][i];
578                b[k]+=v[k]*value;
579              }
580              c += value*value;
581              for (k=0;k<nsolve;k++) {
582                for (j=k;j<nsolve;j++) {
583                  a[k][j] += v[k]*v[j];
584                }
585              }
586            }
587            sum2=c;
588            if (itry==1) {
589              lastObj=objvalue+weight*sum2;
590            }
591            for (k=0;k<nsolve;k++) {
592              b[k] = - (weight*b[k] + 0.5*dj[k]);
593              for (j=k;j<nsolve;j++) {
594                a[k][j] *= weight;
595                a[j][k] = a[k][j];
596              }
597            }
598            c*=weight;
599            for (k=0;k<nsolve;k++) {
600              theta[k]=b[k];
601            }
602            solveSmall(nsolve,a,awork,theta);
603            if ((strategy&64)!=0) {
604              value=10.0;
605              for (k=0;k<nsolve;k++) {
606                value  = max (value, fabs(theta[k]));
607              }
608              if (value>10.0&&((logLevel_&4)!=0)) {
609                printf("theta %g %g %g\n",theta[0],theta[1],theta[2]);
610              }
611              value = 10.0/value;
612              for (k=0;k<nsolve;k++) {
613                theta[k] *= value;
614              }
615            }
616            for (i=0;i<ncolx;i++) {
617              double valueh=0.0;
618              for (k=0;k<HISTORY-1;k++) {
619                value=history[k][i]-history[k+1][i];
620                valueh+=value*theta[k];
621              }
622#ifdef NULLVECTOR
623              value=history[HISTORY][i];
624              valueh+=value*theta[HISTORY-1];
625#endif
626              history[HISTORY][i]=valueh;
627            }
628          }
629#ifdef NULLVECTOR
630          if ((strategy&64)) {
631            for (i=0;i<ncols-nrows;i++) {
632              if (colsol[i]<=lower[i]+tolerance
633                  ||colsol[i]>=(upper[i]-tolerance)) {
634                history[HISTORY][i]=0.0;;
635              }
636            }
637            tolerance=-tolerance; /* switch off test */
638          }
639#endif
640          if (!doScale) {
641            for (i=0;i<ncols;i++) {
642              if (colsol[i]>lower[i]+tolerance
643                  &&colsol[i]<(upper[i]-tolerance)) {
644                value=history[HISTORY][i];
645                colsol[i] +=value;
646                if (colsol[i]<lower[i]+tolerance) {
647                  colsol[i]=lower[i];
648                } else if (colsol[i]>upper[i]-tolerance) {
649                  colsol[i]=upper[i];
650                }
651              }
652            }
653            if (extraBlock) {
654              for (i=0;i<extraBlock;i++) {
655                if (solExtra[i]>tolerance
656                    &&solExtra[i]<(upperExtra[i]-tolerance)) {
657                  value=history[HISTORY][i+ncols];
658                  solExtra[i] +=value;
659                  if (solExtra[i]<0.0) {
660                    solExtra[i]=0.0;
661                  } else if (solExtra[i]>upperExtra[i]) {
662                    solExtra[i]=upperExtra[i];
663                  }
664                }
665              }
666            }
667          } else {
668            double theta=1.0;
669            double saveTheta=theta;
670            for (i=0;i<ncols;i++) {
671              if (colsol[i]>lower[i]+tolerance
672                  &&colsol[i]<(upper[i]-tolerance)) {
673                value=history[HISTORY][i];
674                if (value>0) {
675                  if (theta*value+colsol[i]>upper[i]) {
676                    theta=(upper[i]-colsol[i])/value;
677                  }
678                } else if (value<0) {
679                  if (colsol[i]+theta*value<lower[i]) {
680                    theta = (lower[i]-colsol[i])/value;
681                  }
682                }
683              }
684            }
685            if (extraBlock) {
686              for (i=0;i<extraBlock;i++) {
687                if (solExtra[i]>tolerance
688                    &&solExtra[i]<(upperExtra[i]-tolerance)) {
689                  value=history[HISTORY][i+ncols];
690                  if (value>0) {
691                    if (theta*value+solExtra[i]>upperExtra[i]) {
692                      theta=(upperExtra[i]-solExtra[i])/value;
693                    }
694                  } else if (value<0) {
695                    if (solExtra[i]+theta*value<0.0) {
696                      theta = -solExtra[i]/value;
697                    }
698                  }
699                }
700              }
701            }
702            if ((iter%100==0)&&(logLevel_&8)!=0) {
703              if (theta<saveTheta) {
704                printf(" - modified theta %g\n",theta);
705              }
706            }
707            for (i=0;i<ncols;i++) {
708              if (colsol[i]>lower[i]+tolerance
709                  &&colsol[i]<(upper[i]-tolerance)) {
710                value=history[HISTORY][i];
711                colsol[i] +=value*theta;
712              }
713            }
714            if (extraBlock) {
715              for (i=0;i<extraBlock;i++) {
716                if (solExtra[i]>tolerance
717                    &&solExtra[i]<(upperExtra[i]-tolerance)) {
718                  value=history[HISTORY][i+ncols];
719                  solExtra[i] +=value*theta;
720                }
721              }
722            }
723          }
724#ifdef NULLVECTOR
725          tolerance=fabs(tolerance); /* switch back on */
726#endif
727          if ((iter %100) ==0&&(logLevel_&8)!=0) {
728            printf("\n");
729          }
730        }
731        good=1;
732        result = objval(nrows,ncols,rowsol,colsol,pi,djs,useCost,
733                            rowlower,rowupper,lower,upper,
734                            elemnt,row,columnStart,length,extraBlock,rowExtra,
735                            solExtra,elemExtra,upperExtra,useCostExtra,
736                            weight);
737        weightedObj=result.weighted;
738        if (!iter) saveValue=weightedObj;
739        objvalue=result.objval;
740        sum1=result.infeas;
741        if (saveSol) {
742          if (result.weighted<bestSol) {
743            printf("%d %g better than %g\n",iter,
744                   result.weighted*maxmin-useOffset,bestSol*maxmin-useOffset);
745            bestSol=result.weighted;
746            memcpy(saveSol,colsol,ncols*sizeof(double));
747          }
748        }
749#ifdef FITz
750        if (iter>after) {
751          IdiotResult result2;
752          double ww,oo,ss;
753          if (extraBlock) abort();
754          result2= objval(nrows,ncols,row2,sol2,pi2,djs,cost,
755                     rowlower,rowupper,lower,upper,
756                     elemnt,row,columnStart,extraBlock,rowExtra,
757                     solExtra,elemExtra,upperExtra,costExtra,
758                     weight);
759          ww=result2.weighted;
760          oo=result2.objval;
761          ss=result2.infeas;
762          printf("wobj %g obj %g inf %g last %g\n",ww,oo,ss,lastObj);
763          if (ww<weightedObj&&ww<lastObj) {
764            printf(" taken");
765            ntaken++;
766            saving+=weightedObj-ww;
767            weightedObj=ww;
768            objvalue=oo;
769            sum1=ss;
770            memcpy(rowsol,row2,nrows*sizeof(double));
771            memcpy(pi,pi2,nrows*sizeof(double));
772            memcpy(colsol,sol2,ncols*sizeof(double));
773            result= objval(nrows,ncols,rowsol,colsol,pi,djs,cost,
774                            rowlower,rowupper,lower,upper,
775                            elemnt,row,columnStart,extraBlock,rowExtra,
776                            solExtra,elemExtra,upperExtra,costExtra,
777                            weight);
778            weightedObj=result.weighted;
779            objvalue=result.objval;
780            sum1=result.infeas;
781            if (ww<weightedObj) abort();
782          } else {
783            printf(" not taken");
784            nottaken++;
785          }
786        }
787#endif
788        /*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/
789        if (weightedObj>lastObj+1.0e-4&&itry<MAXTRY) {
790          if((logLevel_&16)!=0&&doScale) {
791            printf("Weighted objective from %g to %g **** bad move\n",
792                   lastObj,weightedObj);
793          }
794          if (doScale) {
795            good=1;
796          }
797          if ((strategy&3)==1) {
798            good=0;
799            if (weightedObj>lastObj+djExit) {
800            if ((logLevel_&16)!=0) {
801              printf("Weighted objective from %g to %g ?\n",lastObj,weightedObj);
802            }
803              memcpy(colsol,history[0],ncols*sizeof(double));
804              memcpy(solExtra,history[0]+ncols,extraBlock*sizeof(double));
805              good=1;
806            }
807          } else if ((strategy&3)==2) {
808            if (weightedObj>lastObj+0.1*maxDj) {
809              memcpy(colsol,history[0],ncols*sizeof(double));
810              memcpy(solExtra,history[0]+ncols,extraBlock*sizeof(double));
811              doScale++;
812              good=0;
813            }
814          } else if ((strategy&3)==3) {
815            if (weightedObj>lastObj+0.001*maxDj) {
816              /*doScale++;*/
817              good=0;
818            }
819          }
820        }
821      }
822      if ((iter % checkFrequency_) ==0) {
823        double best=weightedObj;
824        double test=obj[0];
825        for (i=1;i<DROP;i++) {
826          obj[i-1]=obj[i];
827          if (best>obj[i]) best=obj[i];
828        }
829        obj[DROP-1]=best;
830        if (test-best<drop&&(strategy&8)==0) {
831          if ((logLevel_&8)!=0) {
832          printf("Exiting as drop in %d its is %g after %d iterations\n",
833                 DROP*checkFrequency_,test-best,iter);
834          }
835          goto RETURN;
836        }
837      }
838      if ((iter % logFreq_)==0) {
839        double piSum=0.0;
840        for (i=0;i<nrows;i++) {
841          piSum+=(rowsol[i]+rowupper[i])*pi[i];
842        }
843        if ((logLevel_&2)!=0) {
844        printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n",
845               iter,sum1,objvalue*maxmin-useOffset,weightedObj-useOffset,
846               piSum*maxmin-useOffset,maxDj);
847        }
848      }
849      memcpy(statusWork,statusSave,ncols*sizeof(char));
850      nflagged=0;
851    }
852    nChange=0;
853    doFull=0;
854    maxDj=0.0;
855    // go through forwards or backwards and starting at odd places
856    int itry;
857    for (itry=0;itry<2;itry++) {
858      int icol=start[itry];
859      int istop= stop[itry];
860      for (;icol!=istop;icol += direction) {
861        if (!statusWork[icol]) {
862          CoinBigIndex j;
863          double value=colsol[icol];
864          double djval=cost[icol];
865          double djval2,value2;
866          double theta,a,b,c;
867          if (elemnt) {
868            for (j=columnStart[icol];j<columnStart[icol]+length[icol];j++) {
869              int irow=row[j];
870              djval -= elemnt[j]*pi[irow];
871            }
872          } else {
873            for (j=columnStart[icol];j<columnStart[icol]+length[icol];j++) {
874              int irow=row[j];
875              djval -= pi[irow];
876            }
877          }
878          /*printf("xx iter %d seq %d djval %g value %g\n",
879            iter,i,djval,value);*/
880          if (djval>1.0e-5) {
881            value2=(lower[icol]-value);
882          } else {
883            value2=(upper[icol]-value);
884          }
885          djval2 = djval*value2;
886          djval=fabs(djval);
887          if (djval>djTol) {
888            if (djval2<-1.0e-4) {
889              double valuep,thetap;
890              nChange++;
891              if (djval>maxDj) maxDj=djval;
892              /*if (djval>3.55e6) {
893                printf("big\n");
894                }*/
895              a=0.0;
896              b=0.0;
897              c=0.0;     
898              djval2 = cost[icol];
899              if (elemnt) {
900                for (j=columnStart[icol];j<columnStart[icol]+length[icol];j++) {
901                  int irow=row[j];
902                  double value=rowsol[irow];
903                  c += value*value;
904                  a += elemnt[j]*elemnt[j];
905                  b += value*elemnt[j];
906                }
907              } else {
908                for (j=columnStart[icol];j<columnStart[icol]+length[icol];j++) {
909                  int irow=row[j];
910                  double value=rowsol[irow];
911                  c += value*value;
912                  a += 1.0;
913                  b += value;
914                }
915              }
916              a*=weight;
917              b = b*weight+0.5*djval2;
918              c*=weight;
919              /* solve */
920              theta = -b/a;
921              if ((strategy&4)!=0) {
922                value2=a*theta*theta+2.0*b*theta;
923                thetap=2.0*theta;
924                valuep=a*thetap*thetap+2.0*b*thetap;
925                if (valuep<value2+djTol) {
926                  theta=thetap;
927                  kgood++;
928                } else {
929                  kbad++;
930                }
931              }
932              if (theta>0.0) {
933                if (theta<upper[icol]-colsol[icol]) {
934                  value2=theta;
935                } else {
936                  value2=upper[icol]-colsol[icol];
937                }
938              } else {
939                if (theta>lower[icol]-colsol[icol]) {
940                  value2=theta;
941                } else {
942                  value2=lower[icol]-colsol[icol];
943                }
944              }
945              colsol[icol] += value2;
946              objvalue += cost[icol]*value2;
947              if (elemnt) {
948                for (j=columnStart[icol];j<columnStart[icol]+length[icol];j++) {
949                  int irow=row[j];
950                  double value;
951                  rowsol[irow]+=elemnt[j]*value2;
952                  value=rowsol[irow];
953                  pi[irow]=-2.0*weight*value;
954                }
955              } else {
956                for (j=columnStart[icol];j<columnStart[icol]+length[icol];j++) {
957                  int irow=row[j];
958                  double value;
959                  rowsol[irow]+=value2;
960                  value=rowsol[irow];
961                  pi[irow]=-2.0*weight*value;
962                }
963              }
964            } else {
965              /* dj but at bound */
966              if (djval>djFlag) {
967                statusWork[icol]=1;
968                nflagged++;
969              }
970            }
971          }
972        }
973      }
974    }
975    if (extraBlock) {
976      for (i=0;i<extraBlock;i++) {
977        double value=solExtra[i];
978        double djval=costExtra[i];
979        double djval2,value2;
980        double element=elemExtra[i];
981        double theta,a,b,c;
982        int irow=rowExtra[i];
983        djval -= element*pi[irow];
984        /*printf("xxx iter %d extra %d djval %g value %g\n",
985          iter,irow,djval,value);*/
986        if (djval>1.0e-5) {
987          value2=-value;
988        } else {
989          value2=(upperExtra[i]-value);
990        }
991        djval2 = djval*value2;
992        if (djval2<-1.0e-4&&fabs(djval)>djTol) {
993          nChange++;
994          a=0.0;
995          b=0.0;
996          c=0.0;     
997          djval2 = costExtra[i];
998          value=rowsol[irow];
999          c += value*value;
1000          a += element*element;
1001          b += element*value;
1002          a*=weight;
1003          b = b*weight+0.5*djval2;
1004          c*=weight;
1005          /* solve */
1006          theta = -b/a;
1007          if (theta>0.0) {
1008            value2=min(theta,upperExtra[i]-solExtra[i]);
1009          } else {
1010            value2=max(theta,-solExtra[i]);
1011          }
1012          solExtra[i] += value2;
1013          rowsol[irow]+=element*value2;
1014          value=rowsol[irow];
1015          pi[irow]=-2.0*weight*value;
1016        }
1017      }
1018    }
1019    if ((iter%10)==2) {
1020      for (i=DJTEST-1;i>0;i--) {
1021        djSave[i]=djSave[i-1];
1022      }
1023      djSave[0]=maxDj;
1024      largestDj=max(largestDj,maxDj);
1025      smallestDj=min(smallestDj,maxDj);
1026      for (i=DJTEST-1;i>0;i--) {
1027        maxDj+=djSave[i];
1028      }
1029      maxDj = maxDj/(double) DJTEST;
1030      if (maxDj<djExit&&iter>50) {
1031        //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
1032        break;
1033      }
1034      if (nChange<100) {
1035        djTol*=0.5;
1036      }
1037    }
1038  }
1039 RETURN:
1040  if (kgood||kbad) {
1041    printf("%g good %g bad\n",kgood,kbad);
1042  }
1043  result = objval(nrows,ncols,rowsol,colsol,pi,djs,useCost,
1044                            rowlower,rowupper,lower,upper,
1045                            elemnt,row,columnStart,length,extraBlock,rowExtra,
1046                            solExtra,elemExtra,upperExtra,useCostExtra,
1047                            weight);
1048  result.djAtBeginning=largestDj;
1049  result.djAtEnd=smallestDj;
1050  result.dropThis=saveValue-result.weighted;
1051  if (saveSol) {
1052    if (result.weighted<bestSol) {
1053      bestSol=result.weighted;
1054      memcpy(saveSol,colsol,ncols*sizeof(double));
1055    } else {
1056      printf("restoring previous - now %g best %g\n",
1057             result.weighted*maxmin-useOffset,bestSol*maxmin-useOffset);
1058    }
1059  }
1060  if (saveSol) {
1061    if (extraBlock) {
1062      delete [] useCostExtra;
1063    }
1064    memcpy(colsol,saveSol,ncols*sizeof(double));
1065    delete [] saveSol;
1066  }
1067  for (i=0;i<nsolve;i++) {
1068    if (i) delete [] allsum[i];
1069    delete [] aX[i];
1070    delete [] aworkX[i];
1071  }
1072  delete [] thetaX;
1073  delete [] djX;
1074  delete [] bX;
1075  delete [] vX;
1076  delete [] aX;
1077  delete [] aworkX;
1078  delete [] allsum;
1079  delete [] cost;
1080  for (i=0;i<HISTORY+1;i++) {
1081    delete [] history[i];
1082  }
1083  delete [] statusSave;
1084  /* do original costs objvalue*/
1085  result.objval=0.0;
1086  for (i=0;i<ncols;i++) {
1087    result.objval+=colsol[i]*origcost[i];
1088  }
1089  if (extraBlock) {
1090    for (i=0;i<extraBlock;i++) {
1091      int irow=rowExtra[i];
1092      rowupper[irow]=saveExtra[i];
1093    }
1094    delete [] rowExtra;
1095    delete [] solExtra;
1096    delete [] elemExtra;
1097    delete [] upperExtra;
1098    delete [] costExtra;
1099    delete [] saveExtra;
1100  }
1101  result.iteration=iter;
1102  result.objval-=saveOffset;
1103  result.weighted=result.objval+weight*result.sumSquared;
1104  return result;
1105}
Note: See TracBrowser for help on using the repository browser.