source: branches/devel/Cbc/src/CbcLinked.cpp @ 501

Last change on this file since 501 was 501, checked in by forrest, 14 years ago

qudaratic stuff

File size: 128.0 KB
Line 
1// Copyright (C) 2006, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#include "CbcConfig.h"
4
5#ifdef COIN_HAS_ASL
6#define COIN_HAS_LINK
7#endif
8#ifndef COIN_DEVELOP
9#undef COIN_HAS_LINK
10#endif
11#ifdef COIN_HAS_LINK
12#include <cassert>
13#if defined(_MSC_VER)
14// Turn off compiler warning about long names
15#  pragma warning(disable:4786)
16#endif
17#include "CbcLinked.hpp"
18#include "CoinTime.hpp"
19
20#include "CoinHelperFunctions.hpp"
21#include "CoinIndexedVector.hpp"
22#include "CoinMpsIO.hpp"
23#include "CoinModel.hpp"
24#include "ClpSimplex.hpp"
25//#include "OsiSolverLink.hpp"
26//#include "OsiBranchLink.hpp"
27#include "ClpPackedMatrix.hpp"
28#include "CoinTime.hpp"
29#include "ClpQuadraticObjective.hpp"
30#include "CbcModel.hpp"
31#include "CbcCutGenerator.hpp"
32#include "CglStored.hpp"
33//#include "CglTemporary.hpp"
34//#############################################################################
35// Solve methods
36//#############################################################################
37void OsiSolverLink::initialSolve()
38{
39  //writeMps("yy");
40  //exit(77);
41  specialOptions_ =0;
42  modelPtr_->setWhatsChanged(0);
43  if (numberVariables_) {
44    CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
45    // update all bounds before coefficients
46    for (int i=0;i<numberVariables_;i++ ) {
47      info_[i].updateBounds(modelPtr_);
48    }
49    int updated = updateCoefficients(modelPtr_,temp);
50    if (updated||1) {
51      temp->removeGaps(1.0e-14);
52      ClpMatrixBase * save = modelPtr_->clpMatrix();
53      ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
54      assert (clpMatrix);
55      if (save->getNumRows()>temp->getNumRows()) {
56        // add in cuts
57        int numberRows = temp->getNumRows();
58        int * which = new int[numberRows];
59        for (int i=0;i<numberRows;i++)
60          which[i]=i;
61        save->deleteRows(numberRows,which);
62        delete [] which;
63        temp->bottomAppendPackedMatrix(*clpMatrix->matrix());
64      }
65      modelPtr_->replaceMatrix(temp,true);
66    } else {
67      delete temp;
68    }
69  }
70  if (0) {
71    const double * lower = getColLower();
72    const double * upper = getColUpper();
73    int n=0;
74    for (int i=84;i<84+16;i++) {
75      if (lower[i]+0.01<upper[i]) {
76        n++;
77      }
78    }
79    if (!n)
80      writeMps("sol_query"); 
81
82  }
83  //static int iPass=0;
84  //char temp[50];
85  //iPass++;
86  //sprintf(temp,"cc%d",iPass);
87  //writeMps(temp);
88  //printf("wrote cc%d\n",iPass);
89  OsiClpSolverInterface::initialSolve();
90  int secondaryStatus = modelPtr_->secondaryStatus();
91  if (modelPtr_->status()==0&&(secondaryStatus==2||secondaryStatus==4))
92    modelPtr_->cleanup(1);
93  if (!isProvenOptimal())
94    writeMps("yy");
95  if (isProvenOptimal()&&quadraticModel_&&modelPtr_->numberColumns()==quadraticModel_->numberColumns()) {
96    // see if qp can get better solution
97    const double * solution = modelPtr_->primalColumnSolution();
98    int numberColumns = modelPtr_->numberColumns();
99    bool satisfied=true;
100    for (int i=0;i<numberColumns;i++) {
101      if (isInteger(i)) {
102        double value = solution[i];
103        if (fabs(value-floor(value+0.5))>1.0e-6) {
104          satisfied=false;
105          break;
106        }
107      }
108    }
109    if (satisfied) {
110      ClpSimplex qpTemp(*quadraticModel_);
111      double * lower = qpTemp.columnLower();
112      double * upper = qpTemp.columnUpper();
113      double * lower2 = modelPtr_->columnLower();
114      double * upper2 = modelPtr_->columnUpper();
115      for (int i=0;i<numberColumns;i++) {
116        if (isInteger(i)) {
117          double value = floor(solution[i]+0.5);
118          lower[i]=value;
119          upper[i]=value;
120        } else {
121          lower[i]=lower2[i];
122          upper[i]=upper2[i];
123        }
124      }
125      //qpTemp.writeMps("bad.mps");
126      //modelPtr_->writeMps("bad2.mps");
127      //qpTemp.objectiveAsObject()->setActivated(0);
128      //qpTemp.primal();
129      //qpTemp.objectiveAsObject()->setActivated(1);
130      qpTemp.primal();
131      //assert (!qpTemp.problemStatus());
132      if (qpTemp.objectiveValue()<bestObjectiveValue_-1.0e-3&&!qpTemp.problemStatus()) {
133        delete [] bestSolution_;
134        bestSolution_ = CoinCopyOfArray(qpTemp.primalColumnSolution(),numberColumns);
135        bestObjectiveValue_ = qpTemp.objectiveValue();
136        printf("better qp objective of %g\n",bestObjectiveValue_);
137        // If model has stored then add cut (if convex)
138        if (cbcModel_&&(specialOptions2_&4)!=0) {
139          int numberGenerators = cbcModel_->numberCutGenerators();
140          int iGenerator;
141          for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
142            CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
143            CglCutGenerator * gen = generator->generator();
144            CglStored * gen2 = dynamic_cast<CglStored *> (gen);
145            if (gen2) {
146              // add OA cut
147              double offset;
148              double * gradient = new double [numberColumns+1];
149              memcpy(gradient,qpTemp.objectiveAsObject()->gradient(&qpTemp,bestSolution_,offset,true,2),
150                     numberColumns*sizeof(double));
151              // assume convex
152              double rhs = 0.0;
153              int * column = new int[numberColumns+1];
154              int n=0;
155              for (int i=0;i<numberColumns;i++) {
156                double value = gradient[i];
157                if (fabs(value)>1.0e-12) {
158                  gradient[n]=value;
159                  rhs += value*solution[i];
160                  column[n++]=i;
161                }
162              }
163              gradient[n]=-1.0;
164              column[n++]=objectiveVariable_;
165              gen2->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
166              delete [] gradient;
167              delete [] column;
168              break;
169            }
170          }
171        }
172      }
173    }
174  }
175}
176//#define WRITE_MATRIX
177#ifdef WRITE_MATRIX
178static int xxxxxx=0;
179#endif
180//-----------------------------------------------------------------------------
181void OsiSolverLink::resolve()
182{
183  specialOptions_ =0;
184  modelPtr_->setWhatsChanged(0);
185  bool allFixed=false;
186  bool feasible=true;
187  if (numberVariables_) {
188    CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
189    allFixed=true;
190    //bool best=true;
191    const double * lower = modelPtr_->columnLower();
192    const double * upper = modelPtr_->columnUpper();
193    // update all bounds before coefficients
194    for (int i=0;i<numberVariables_;i++ ) {
195      info_[i].updateBounds(modelPtr_);
196      int iColumn = info_[i].variable();
197      double lo = lower[iColumn];
198      double up = upper[iColumn];
199      if (up>lo)
200        allFixed=false;
201      else if (up<lo)
202        feasible=false;
203    }
204    int updated=updateCoefficients(modelPtr_,temp);
205    if (updated) {
206      temp->removeGaps(1.0e-14);
207      ClpMatrixBase * save = modelPtr_->clpMatrix();
208      ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
209      assert (clpMatrix);
210      if (save->getNumRows()>temp->getNumRows()) {
211        // add in cuts
212        int numberRows = temp->getNumRows();
213        int * which = new int[numberRows];
214        for (int i=0;i<numberRows;i++)
215          which[i]=i;
216        CoinPackedMatrix * mat = clpMatrix->matrix();
217        // for debug
218        //mat = new CoinPackedMatrix(*mat);
219        mat->deleteRows(numberRows,which);
220        delete [] which;
221        temp->bottomAppendPackedMatrix(*mat);
222        temp->removeGaps(1.0e-14);
223      }
224      if (0) {
225        const CoinPackedMatrix * matrix = modelPtr_->matrix();
226        int numberColumns = matrix->getNumCols();
227        assert (numberColumns==temp->getNumCols());
228        const double * element1 = temp->getMutableElements();
229        const int * row1 = temp->getIndices();
230        const CoinBigIndex * columnStart1 = temp->getVectorStarts();
231        const int * columnLength1 = temp->getVectorLengths();
232        const double * element2 = matrix->getMutableElements();
233        const int * row2 = matrix->getIndices();
234        const CoinBigIndex * columnStart2 = matrix->getVectorStarts();
235        const int * columnLength2 = matrix->getVectorLengths();
236        for (int i=0;i<numberColumns;i++) {
237          assert (columnLength2[i]==columnLength1[i]);
238          int offset = columnStart2[i]-columnStart1[i];
239          for (int j=columnStart1[i];j<columnStart1[i]+columnLength1[i];j++) {
240            assert (row1[j]==row2[j+offset]);
241            assert (element1[j]==element2[j+offset]);
242          }
243        }
244      }
245      modelPtr_->replaceMatrix(temp,true);
246    } else {
247      delete temp;
248    }
249  }
250#ifdef WRITE_MATRIX
251  {
252    xxxxxx++;
253    char temp[50];
254    sprintf(temp,"bb%d",xxxxxx);
255    writeMps(temp);
256    printf("wrote bb%d\n",xxxxxx);
257#if 0
258    const double * lower = getColLower();
259    const double * upper = getColUpper();
260    for (int i=0;i<8;i++) {
261      printf("%d bounds %g %g\n",i,lower[i],upper[i]);
262    }
263    if (lower[2]<=3.0&&upper[2]>=3) {
264      if (lower[3]<=2.0&&upper[3]>=2.0) {
265        if (lower[4]==0.0) {
266          if (lower[7]==0.0) {
267            if (lower[5]<=4.0&&upper[5]>=4.0) {
268              if (lower[6]<=3.0&&upper[6]>=3.0) {
269                printf("writemps contains solution\n");
270                for (int i=8;i<24;i++) {
271                  printf("%d bounds %g %g\n",i,lower[i],upper[i]);
272                  if (false&&((!lower[i]&&upper[i]>lower[i]))) {
273                    setColLower(i,0.0);
274                    setColUpper(i,10.0);
275                  }
276                }
277              }
278            }
279          }
280        }
281      }
282    }
283#endif
284  }
285#endif
286  if (0) {
287    const double * lower = getColLower();
288    const double * upper = getColUpper();
289    int n=0;
290    for (int i=60;i<64;i++) {
291      if (lower[i]) {
292        printf("%d bounds %g %g\n",i,lower[i],upper[i]);
293        n++;
294      }
295    }
296    if (n==1)
297      printf("just one?\n");
298  }
299  // check feasible
300   {
301    const double * lower = getColLower();
302    const double * upper = getColUpper();
303    int numberColumns = getNumCols();
304    for (int i=0;i<numberColumns;i++) {
305      if (lower[i]>upper[i]+1.0e-12) {
306        feasible = false;
307        break;
308      }
309    }
310  }
311  if (!feasible)
312    allFixed=false;
313  if ((specialOptions2_&1)!=0)
314    allFixed=false;
315  int returnCode=-1;
316  // See if in strong branching
317  int maxIts = modelPtr_->maximumIterations();
318  if (feasible) {
319    if(maxIts>10000) {
320      // may do lots of work
321      returnCode=fathom(allFixed);
322    } else {
323      returnCode=0;
324    }
325  }
326  if (returnCode>=0) {
327    if (returnCode==0)
328      OsiClpSolverInterface::resolve();
329    if (!allFixed&&(specialOptions2_&1)==0) {
330      const double * solution = getColSolution();
331      bool satisfied=true;
332      for (int i=0;i<numberVariables_;i++) {
333        int iColumn = info_[i].variable();
334        double value = solution[iColumn];
335        if (fabs(value-floor(value+0.5))>0.0001)
336          satisfied=false;
337      }
338      //if (satisfied)
339      //printf("satisfied but not fixed\n");
340    }
341    int satisfied=2;
342    const double * solution = getColSolution();
343    const double * lower = getColLower();
344    const double * upper = getColUpper();
345    int numberColumns2 = coinModel_.numberColumns();
346    for (int i=0;i<numberColumns2;i++) {
347      if (isInteger(i)) {
348        double value = solution[i];
349        if (fabs(value-floor(value+0.5))>1.0e-6) {
350          satisfied=0;
351          break;
352        } else if (upper[i]>lower[i]) {
353          satisfied=1;
354        }
355      }
356    }
357    if (isProvenOptimal()) {
358      if (satisfied&&(specialOptions2_&2)!=0) {
359        assert (quadraticModel_);
360        // look at true objective
361        double direction = modelPtr_->optimizationDirection();
362        assert (direction==1.0);
363        double value = - quadraticModel_->objectiveOffset();
364        const double * objective = quadraticModel_->objective();
365        int i;
366        for ( i=0;i<numberColumns2;i++) 
367          value += solution[i]*objective[i];
368        // and now rest
369        for ( i =0;i<numberObjects_;i++) {
370          OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
371          if (obj) {
372            value += obj->xyCoefficient(solution);
373          }
374        }
375        if (value<bestObjectiveValue_-1.0e-3) {
376          printf("obj of %g\n",value);
377          //modelPtr_->setDualObjectiveLimit(value);
378          delete [] bestSolution_;
379          bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(),modelPtr_->getNumCols());
380          bestObjectiveValue_ = value;
381          // If model has stored then add cut (if convex)
382          if (cbcModel_&&(specialOptions2_&4)!=0&&quadraticModel_) {
383            int numberGenerators = cbcModel_->numberCutGenerators();
384            int iGenerator;
385            for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
386              CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
387              CglCutGenerator * gen = generator->generator();
388              CglStored * gen2 = dynamic_cast<CglStored *> (gen);
389              if (gen2) {
390                // add OA cut
391                double offset;
392                int numberColumns = quadraticModel_->numberColumns();
393                double * gradient = new double [numberColumns+1];
394                // gradient from bilinear
395                int i;
396                CoinZeroN(gradient,numberColumns+1);
397                //const double * objective = modelPtr_->objective();
398                assert (objectiveRow_>=0);
399                const double * element = originalRowCopy_->getElements();
400                const int * column2 = originalRowCopy_->getIndices();
401                const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
402                //const int * rowLength = originalRowCopy_->getVectorLengths();
403                //int numberColumns2 = coinModel_.numberColumns();
404                for ( i=rowStart[objectiveRow_];i<rowStart[objectiveRow_+1];i++) 
405                  gradient[column2[i]] = element[i];
406                for ( i =0;i<numberObjects_;i++) {
407                  OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
408                  if (obj) {
409                    int xColumn = obj->xColumn();
410                    int yColumn = obj->yColumn();
411                    if (xColumn!=yColumn) {
412                      double coefficient = 2.0*obj->coefficient();
413                      gradient[xColumn] += coefficient*solution[yColumn];
414                      gradient[yColumn] += coefficient*solution[xColumn];
415                      offset += coefficient*solution[xColumn]*solution[yColumn];
416                    } else {
417                      double coefficient = obj->coefficient();
418                      gradient[xColumn] += 2.0*coefficient*solution[yColumn];
419                      offset += coefficient*solution[xColumn]*solution[yColumn];
420                    }
421                  }
422                }
423                // assume convex
424                double rhs = 0.0;
425                int * column = new int[numberColumns+1];
426                int n=0;
427                for (int i=0;i<numberColumns;i++) {
428                  double value = gradient[i];
429                  if (fabs(value)>1.0e-12) {
430                    gradient[n]=value;
431                    rhs += value*solution[i];
432                    column[n++]=i;
433                  }
434                }
435                gradient[n]=-1.0;
436                column[n++]=objectiveVariable_;
437                gen2->addCut(-COIN_DBL_MAX,offset+1.0e-4,n,column,gradient);
438                delete [] gradient;
439                delete [] column;
440                break;
441              }
442            }
443          }
444        }
445      } else if (satisfied==2) {
446        // is there anything left to do?
447        int i;
448        int numberContinuous=0;
449        double gap=0.0;
450        for ( i =0;i<numberObjects_;i++) {
451          OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
452          if (obj) {
453            if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
454              numberContinuous++;
455              int xColumn = obj->xColumn();
456              double gapX = upper[xColumn]-lower[xColumn];
457              int yColumn = obj->yColumn();
458              double gapY = upper[yColumn]-lower[yColumn];
459              gap = CoinMax(gap,CoinMax(gapX,gapY));
460            }
461          }
462        }
463        if (numberContinuous&&0) {
464          // iterate to get solution and fathom node
465          int numberColumns2 = coinModel_.numberColumns();
466          double * lower2 = CoinCopyOfArray(getColLower(),numberColumns2);
467          double * upper2 = CoinCopyOfArray(getColUpper(),numberColumns2);
468          while (gap>defaultMeshSize_) {
469            gap *= 0.9;
470            const double * solution = getColSolution();
471            double newGap=0.0;
472            for ( i =0;i<numberObjects_;i++) {
473              OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
474              if (obj&&(obj->branchingStrategy()&8)==0) {
475                if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0) {
476                  numberContinuous++;
477                  // need to make sure new xy value in range
478                  double xB[3];
479                  double yB[3];
480                  double xybar[4];
481                  obj->getCoefficients(this,xB,yB,xybar);
482                  //double xyTrue = obj->xyCoefficient(solution);
483                  double xyLambda=0.0;
484                  int firstLambda = obj->firstLambda();
485                  for (int j=0;j<4;j++) {
486                    xyLambda += solution[firstLambda+j]*xybar[j];
487                  }
488                  //printf ("x %d, y %d - true %g lambda %g\n",obj->xColumn(),obj->yColumn(),
489                  //  xyTrue,xyLambda);
490                  int xColumn = obj->xColumn();
491                  double gapX = upper[xColumn]-lower[xColumn];
492                  int yColumn = obj->yColumn();
493                  if (gapX>gap) {
494                    double value = solution[xColumn];
495                    double newLower = CoinMax(lower2[xColumn],value-0.5*gap);
496                    double newUpper = CoinMin(upper2[xColumn],value+0.5*gap);
497                    if (newUpper-newLower<0.99*gap) {
498                      if (newLower==lower2[xColumn])
499                        newUpper = CoinMin(upper2[xColumn],newLower+gap);
500                      else if (newUpper==upper2[xColumn])
501                        newLower = CoinMax(lower2[xColumn],newUpper-gap);
502                    }
503                    // see if problem
504                    double lambda[4];
505                    xB[0]=newLower;
506                    xB[1]=newUpper;
507                    xB[2]=value;
508                    yB[2]=solution[yColumn];
509                    xybar[0]=xB[0]*yB[0];
510                    xybar[1]=xB[0]*yB[1];
511                    xybar[2]=xB[1]*yB[0];
512                    xybar[3]=xB[1]*yB[1];
513                    double infeasibility=obj->computeLambdas(xB,yB,xybar,lambda);
514                    assert (infeasibility<1.0e-9);
515                    setColLower(xColumn,newLower);
516                    setColUpper(xColumn,newUpper);
517                  }
518                  double gapY = upper[yColumn]-lower[yColumn];
519                  if (gapY>gap) {
520                    double value = solution[yColumn];
521                    double newLower = CoinMax(lower2[yColumn],value-0.5*gap);
522                    double newUpper = CoinMin(upper2[yColumn],value+0.5*gap);
523                    if (newUpper-newLower<0.99*gap) {
524                      if (newLower==lower2[yColumn])
525                        newUpper = CoinMin(upper2[yColumn],newLower+gap);
526                      else if (newUpper==upper2[yColumn])
527                        newLower = CoinMax(lower2[yColumn],newUpper-gap);
528                    }
529                    // see if problem
530                    double lambda[4];
531                    yB[0]=newLower;
532                    yB[1]=newUpper;
533                    xybar[0]=xB[0]*yB[0];
534                    xybar[1]=xB[0]*yB[1];
535                    xybar[2]=xB[1]*yB[0];
536                    xybar[3]=xB[1]*yB[1];
537                    double infeasibility=obj->computeLambdas(xB,yB,xybar,lambda);
538                    assert (infeasibility<1.0e-9);
539                    setColLower(yColumn,newLower);
540                    setColUpper(yColumn,newUpper);
541                  }
542                  newGap = CoinMax(newGap,CoinMax(gapX,gapY));
543                }
544              }
545            }
546            printf("solving with gap of %g\n",gap);
547            //OsiClpSolverInterface::resolve();
548            initialSolve();
549            if (!isProvenOptimal()) 
550              break;
551          }
552          delete [] lower2;
553          delete [] upper2;
554          if (isProvenOptimal())
555            writeMps("zz");
556        }
557      }
558      // ???  - try
559      // But skip if strong branching
560      CbcModel * cbcModel = (modelPtr_->maximumIterations()<10000) ? cbcModel_ : NULL;
561      if ((specialOptions2_&2)!=0) {
562        // If model has stored then add cut (if convex)
563        if (cbcModel&&(specialOptions2_&4)!=0&&quadraticModel_) {
564          int numberGenerators = cbcModel_->numberCutGenerators();
565          int iGenerator;
566          for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
567            CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
568            CglCutGenerator * gen = generator->generator();
569            CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
570            if (gen2) {
571              const double * solution = getColSolution();
572              // add OA cut
573              double offset;
574              int numberColumns = quadraticModel_->numberColumns();
575              double * gradient = new double [numberColumns+1];
576              // gradient from bilinear
577              int i;
578              CoinZeroN(gradient,numberColumns+1);
579              //const double * objective = modelPtr_->objective();
580              assert (objectiveRow_>=0);
581              const double * element = originalRowCopy_->getElements();
582              const int * column2 = originalRowCopy_->getIndices();
583              const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
584              //const int * rowLength = originalRowCopy_->getVectorLengths();
585              //int numberColumns2 = coinModel_.numberColumns();
586              for ( i=rowStart[objectiveRow_];i<rowStart[objectiveRow_+1];i++) 
587                gradient[column2[i]] = element[i];
588              for ( i =0;i<numberObjects_;i++) {
589                OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
590                if (obj) {
591                  int xColumn = obj->xColumn();
592                  int yColumn = obj->yColumn();
593                  if (xColumn!=yColumn) {
594                    double coefficient = 2.0*obj->coefficient();
595                    gradient[xColumn] += coefficient*solution[yColumn];
596                    gradient[yColumn] += coefficient*solution[xColumn];
597                    offset += coefficient*solution[xColumn]*solution[yColumn];
598                  } else {
599                    double coefficient = obj->coefficient();
600                    gradient[xColumn] += 2.0*coefficient*solution[yColumn];
601                    offset += coefficient*solution[xColumn]*solution[yColumn];
602                  }
603                }
604              }
605              // assume convex
606              double rhs = 0.0;
607              int * column = new int[numberColumns+1];
608              int n=0;
609              for (int i=0;i<numberColumns;i++) {
610                double value = gradient[i];
611                if (fabs(value)>1.0e-12) {
612                  gradient[n]=value;
613                  rhs += value*solution[i];
614                  column[n++]=i;
615                }
616              }
617              gradient[n]=-1.0;
618              assert (objectiveVariable_>=0);
619              rhs -= solution[objectiveVariable_];
620              column[n++]=objectiveVariable_;
621              if (rhs>offset+1.0e-5) {
622                gen2->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
623                //printf("added cut with %d elements\n",n);
624              }
625              delete [] gradient;
626              delete [] column;
627              break;
628            }
629          }
630        }
631      } else if (cbcModel&&(specialOptions2_&8)==8) {
632        // convex and nonlinear in constraints
633        int numberGenerators = cbcModel_->numberCutGenerators();
634        int iGenerator;
635        for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
636          CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
637          CglCutGenerator * gen = generator->generator();
638          CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
639          if (gen2) {
640            const double * solution = getColSolution();
641            const double * rowUpper = getRowUpper();
642            const double * rowLower = getRowLower();
643            const double * element = originalRowCopy_->getElements();
644            const int * column2 = originalRowCopy_->getIndices();
645            const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
646            //const int * rowLength = originalRowCopy_->getVectorLengths();
647            int numberColumns2 = CoinMax(coinModel_.numberColumns(),objectiveVariable_+1);
648            double * gradient = new double [numberColumns2];
649            int * column = new int[numberColumns2];
650            for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
651              int iRow = rowNonLinear_[iNon];
652              bool convex = convex_[iNon]>0;
653              if (!convex_[iNon])
654                continue; // can't use this row
655              // add OA cuts
656              double offset=0.0;
657              // gradient from bilinear
658              int i;
659              CoinZeroN(gradient,numberColumns2);
660              //const double * objective = modelPtr_->objective();
661              for ( i=rowStart[iRow];i<rowStart[iRow+1];i++) 
662                gradient[column2[i]] = element[i];
663              for ( i =startNonLinear_[iNon];i<startNonLinear_[iNon+1];i++) {
664                OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
665                assert (obj);
666                int xColumn = obj->xColumn();
667                int yColumn = obj->yColumn();
668                if (xColumn!=yColumn) {
669                  double coefficient = 2.0*obj->coefficient();
670                  gradient[xColumn] += coefficient*solution[yColumn];
671                  gradient[yColumn] += coefficient*solution[xColumn];
672                  offset += coefficient*solution[xColumn]*solution[yColumn];
673                } else {
674                  double coefficient = obj->coefficient();
675                  gradient[xColumn] += 2.0*coefficient*solution[yColumn];
676                  offset += coefficient*solution[xColumn]*solution[yColumn];
677                }
678              }
679              // assume convex
680              double rhs = 0.0;
681              int n=0;
682              for (int i=0;i<numberColumns2;i++) {
683                double value = gradient[i];
684                if (fabs(value)>1.0e-12) {
685                  gradient[n]=value;
686                  rhs += value*solution[i];
687                  column[n++]=i;
688                }
689              }
690              if (iRow==objectiveRow_) {
691                gradient[n]=-1.0;
692                assert (objectiveVariable_>=0);
693                rhs -= solution[objectiveVariable_];
694                column[n++]=objectiveVariable_;
695                assert (convex);
696              } else if (convex) {
697                offset += rowUpper[iRow];
698              } else if (!convex) {
699                offset += rowLower[iRow];
700              }
701              if (convex&&rhs>offset+1.0e-5) 
702                gen2->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
703              else if (!convex&&rhs<offset-1.0e-5) 
704                gen2->addCut(offset-1.0e-7,COIN_DBL_MAX,n,column,gradient);
705            }
706            delete [] gradient;
707            delete [] column;
708            break;
709          }
710        }
711      }
712    }
713  } else {
714    modelPtr_->setProblemStatus(1);
715    modelPtr_->setObjectiveValue(COIN_DBL_MAX);
716  }
717}
718
719//#############################################################################
720// Constructors, destructors clone and assignment
721//#############################################################################
722
723//-------------------------------------------------------------------
724// Default Constructor
725//-------------------------------------------------------------------
726OsiSolverLink::OsiSolverLink ()
727  : OsiClpSolverInterface()
728{
729  gutsOfDestructor(true);
730}
731#if 0
732/* returns
733   sequence of nonlinear or
734   -1 numeric
735   -2 not found
736   -3 too many terms
737*/
738static int getVariable(const CoinModel & model, char * expression,
739                       int & linear)
740{
741  int non=-1;
742  linear=-1;
743  if (strcmp(expression,"Numeric")) {
744    // function
745    char * first = strchr(expression,'*');
746    int numberColumns = model.numberColumns();
747    int j;
748    if (first) {
749      *first='\0';
750      for (j=0;j<numberColumns;j++) {
751        if (!strcmp(expression,model.columnName(j))) {
752          linear=j;
753          memmove(expression,first+1,strlen(first+1)+1);
754          break;
755        }
756      }
757    }
758    // find nonlinear
759    for (j=0;j<numberColumns;j++) {
760      const char * name = model.columnName(j);
761      first = strstr(expression,name);
762      if (first) {
763        if (first!=expression&&isalnum(*(first-1)))
764          continue; // not real match
765        first += strlen(name);
766        if (!isalnum(*first)) {
767          // match
768          non=j;
769          // but check no others
770          j++;
771          for (;j<numberColumns;j++) {
772            const char * name = model.columnName(j);
773            first = strstr(expression,name);
774            if (first) {
775              if (isalnum(*(first-1)))
776                continue; // not real match
777              first += strlen(name);
778              if (!isalnum(*first)) {
779                // match - ouch
780                non=-3;
781                break;
782              }
783            }
784          }
785          break;
786        }
787      }
788    }
789    if (non==-1)
790      non=-2;
791  }
792  return non;
793}
794#endif
795/* This creates from a coinModel object
796
797   if errors.then number of sets is -1
798     
799   This creates linked ordered sets information.  It assumes -
800
801   for product terms syntax is yy*f(zz)
802   also just f(zz) is allowed
803   and even a constant
804
805   modelObject not const as may be changed as part of process.
806*/
807OsiSolverLink::OsiSolverLink ( CoinModel & coinModel)
808  : OsiClpSolverInterface()
809{
810  gutsOfDestructor(true);
811  load(coinModel);
812}
813// need bounds
814static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue)
815{
816  double lo = solver->getColLower()[column];
817  if (lo<-maximumValue)
818    solver->setColLower(column,-maximumValue);
819  double up = solver->getColUpper()[column];
820  if (up>maximumValue)
821    solver->setColUpper(column,maximumValue);
822}
823// returns jColumn (-2 if linear term, -1 if unknown) and coefficient
824static
825int decodeBit(char * phrase, char * & nextPhrase, double & coefficient, bool ifFirst, const CoinModel & model)
826{
827  char * pos = phrase;
828  // may be leading - (or +)
829  char * pos2 = pos;
830  double value=1.0;
831  if (*pos2=='-'||*pos2=='+')
832    pos2++;
833  // next terminator * or + or -
834  while (*pos2) {
835    if (*pos2=='*') {
836      break;
837    } else if (*pos2=='-'||*pos2=='+') {
838      if (pos2==pos||*(pos2-1)!='e')
839        break;
840    }
841    pos2++;
842  }
843  // if * must be number otherwise must be name
844  if (*pos2=='*') {
845    char * pos3 = pos;
846    while (pos3!=pos2) {
847      char x = *pos3;
848      pos3++;
849      assert ((x>='0'&&x<='9')||x=='.'||x=='+'||x=='-'||x=='e');
850    }
851    char saved = *pos2;
852    *pos2='\0';
853    value = atof(pos);
854    *pos2=saved;
855    // and down to next
856    pos2++;
857    pos=pos2;
858    while (*pos2) {
859      if (*pos2=='-'||*pos2=='+')
860        break;
861      pos2++;
862    }
863  }
864  char saved = *pos2;
865  *pos2='\0';
866  // now name
867  // might have + or -
868  if (*pos=='+') {
869    pos++;
870  } else if (*pos=='-') {
871    pos++;
872    assert (value==1.0);
873    value = - value;
874  }
875  int jColumn = model.column(pos);
876  // must be column unless first when may be linear term
877  if (jColumn<0) {
878    if (ifFirst) {
879      char * pos3 = pos;
880      while (pos3!=pos2) {
881        char x = *pos3;
882        pos3++;
883        assert ((x>='0'&&x<='9')||x=='.'||x=='+'||x=='-'||x=='e');
884      }
885      assert(*pos2=='\0');
886      // keep possible -
887      value = value * atof(pos);
888      jColumn=-2;
889    } else {
890      // bad
891      *pos2=saved;
892      printf("bad nonlinear term %s\n",phrase);
893      abort();
894    }
895  }
896  *pos2=saved;
897  pos=pos2;
898  coefficient=value;
899  nextPhrase = pos;
900  return jColumn;
901}
902void OsiSolverLink::load ( CoinModel & coinModel, bool tightenBounds)
903{
904  // first check and set up arrays
905  int numberColumns = coinModel.numberColumns();
906  int numberRows = coinModel.numberRows();
907  // List of nonlinear entries
908  int * which = new int[numberColumns];
909  numberVariables_=0;
910  //specialOptions2_=0;
911  int iColumn;
912  int numberErrors=0;
913  for (iColumn=0;iColumn<numberColumns;iColumn++) {
914    CoinModelLink triple=coinModel.firstInColumn(iColumn);
915    bool linear=true;
916    int n=0;
917    // See if quadratic objective
918    const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
919    if (strcmp(expr,"Numeric")) {
920      linear=false;
921    }
922    while (triple.row()>=0) {
923      int iRow = triple.row();
924      const char * expr = coinModel.getElementAsString(iRow,iColumn);
925      if (strcmp(expr,"Numeric")) {
926        linear=false;
927      }
928      triple=coinModel.next(triple);
929      n++;
930    }
931    if (!linear) {
932      which[numberVariables_++]=iColumn;
933    }
934  }
935  // return if nothing
936  if (!numberVariables_) {
937    delete [] which;
938    coinModel_ = coinModel;
939    int nInt=0;
940    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
941      if (coinModel_.isInteger(iColumn))
942        nInt++;
943    }
944    printf("There are %d integers\n",nInt);
945    loadFromCoinModel(coinModel,true);
946    OsiObject ** objects = new OsiObject * [nInt];
947    nInt=0;
948    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
949      if (coinModel_.isInteger(iColumn)) {
950        objects[nInt] = new OsiSimpleInteger(this,iColumn);
951        objects[nInt]->setPriority(integerPriority_);
952        nInt++;
953      }
954    }
955    addObjects(nInt,objects);
956    int i;
957    for (i=0;i<nInt;i++)
958      delete objects[i];
959    delete [] objects;
960    return;
961  } else {
962    coinModel_ = coinModel;
963    // arrays for tightening bounds
964    int * freeRow = new int [numberRows];
965    CoinZeroN(freeRow,numberRows);
966    int * tryColumn = new int [numberColumns];
967    CoinZeroN(tryColumn,numberColumns);
968    int nBi=0;
969    int numberQuadratic=0;
970    for (iColumn=0;iColumn<numberColumns;iColumn++) {
971      // See if quadratic objective
972      const char * expr = coinModel_.getColumnObjectiveAsString(iColumn);
973      if (strcmp(expr,"Numeric")) {
974        // check if value*x+-value*y....
975        assert (strlen(expr)<20000);
976        tryColumn[iColumn]=1;
977        char temp[20000];
978        strcpy(temp,expr);
979        char * pos = temp;
980        bool ifFirst=true;
981        double linearTerm=0.0;
982        while (*pos) {
983          double value;
984          int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
985          // must be column unless first when may be linear term
986          if (jColumn>=0) {
987            tryColumn[jColumn]=1;
988            numberQuadratic++;
989            nBi++;
990          } else if (jColumn==-2) {
991            linearTerm = value;
992          } else {
993            printf("bad nonlinear term %s\n",temp);
994            abort();
995          }
996          ifFirst=false;
997        }
998        coinModel.setObjective(iColumn,linearTerm);
999      }
1000    }
1001    int iRow;
1002    int saveNBi=nBi;
1003    for (iRow=0;iRow<numberRows;iRow++) {   
1004      CoinModelLink triple=coinModel_.firstInRow(iRow);
1005      while (triple.column()>=0) {
1006        int iColumn = triple.column();
1007        const char *  el = coinModel_.getElementAsString(iRow,iColumn);
1008        if (strcmp("Numeric",el)) {
1009          // check if value*x+-value*y....
1010          assert (strlen(el)<20000);
1011          char temp[20000];
1012          strcpy(temp,el);
1013          char * pos = temp;
1014          bool ifFirst=true;
1015          double linearTerm=0.0;
1016          tryColumn[iColumn]=1;
1017          freeRow[iRow]=1;
1018          while (*pos) {
1019            double value;
1020            int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1021            // must be column unless first when may be linear term
1022            if (jColumn>=0) {
1023              tryColumn[jColumn]=1;
1024              nBi++;
1025            } else if (jColumn==-2) {
1026              linearTerm = value;
1027            } else {
1028              printf("bad nonlinear term %s\n",temp);
1029              abort();
1030            }
1031            ifFirst=false;
1032          }
1033          coinModel.setElement(iRow,iColumn,linearTerm);
1034        }
1035        triple=coinModel_.next(triple);
1036      }
1037    }
1038    if (!nBi)
1039      exit(1);
1040    bool quadraticObjective=false;
1041    int nInt=0;
1042    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
1043      if (coinModel_.isInteger(iColumn))
1044        nInt++;
1045    }
1046    printf("There are %d bilinear and %d integers\n",nBi,nInt);
1047    loadFromCoinModel(coinModel,true);
1048    if (tightenBounds) {
1049      // first fake bounds
1050      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1051        if (tryColumn[iColumn]) {
1052          fakeBounds(this,iColumn,defaultBound_);
1053        }
1054      }
1055      ClpSimplex tempModel(*modelPtr_);
1056      int nDelete = 0;
1057      for (iRow=0;iRow<numberRows;iRow++) {
1058        if (freeRow[iRow]) 
1059          freeRow[nDelete++]=iRow;
1060      }
1061      tempModel.deleteRows(nDelete,freeRow);
1062      tempModel.setOptimizationDirection(1.0);
1063      double * objective = tempModel.objective();
1064      CoinZeroN(objective,numberColumns);
1065      // now up and down
1066      double * columnLower = modelPtr_->columnLower();
1067      double * columnUpper = modelPtr_->columnUpper();
1068      const double * solution = tempModel.primalColumnSolution();
1069      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1070        if (tryColumn[iColumn]) {
1071          objective[iColumn]=1.0;
1072          tempModel.primal(1);
1073          if (solution[iColumn]>columnLower[iColumn]+1.0e-3) {
1074            double value = solution[iColumn];
1075            if (coinModel_.isInteger(iColumn))
1076              value = ceil(value-0.9e-3);
1077            printf("lower bound on %d changed from %g to %g\n",iColumn,columnLower[iColumn],value);
1078            columnLower[iColumn]=value;
1079          }
1080          objective[iColumn]=-1.0;
1081          tempModel.primal(1);
1082          if (solution[iColumn]<columnUpper[iColumn]-1.0e-3) {
1083            double value = solution[iColumn];
1084            if (coinModel_.isInteger(iColumn))
1085              value = floor(value+0.9e-3);
1086            printf("upper bound on %d changed from %g to %g\n",iColumn,columnUpper[iColumn],value);
1087            columnUpper[iColumn]=value;
1088          }
1089          objective[iColumn]=0.0;
1090        }
1091      }
1092    }
1093    delete [] freeRow;
1094    delete [] tryColumn;
1095    CoinBigIndex * startQuadratic = NULL;
1096    int * columnQuadratic = NULL;
1097    double * elementQuadratic = NULL;
1098    if( saveNBi==nBi) {
1099      printf("all bilinearity in objective\n");
1100      specialOptions2_ |= 2;
1101      quadraticObjective = true;
1102      // save copy as quadratic model
1103      quadraticModel_ = new ClpSimplex(*modelPtr_);
1104      startQuadratic = new CoinBigIndex [numberColumns+1];
1105      columnQuadratic = new int [numberQuadratic];
1106      elementQuadratic = new double [numberQuadratic];
1107      numberQuadratic=0;
1108    }
1109    if (quadraticObjective||((specialOptions2_&8)!=0&&saveNBi)) {
1110      // add in objective as constraint
1111      objectiveVariable_= numberColumns;
1112      objectiveRow_ = modelPtr_->numberRows();
1113      addCol(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0);
1114      int * column = new int[numberColumns+1];
1115      double * element = new double[numberColumns+1];
1116      double * objective = modelPtr_->objective();
1117      int n=0;
1118      int starts[2]={0,0};
1119      for (int i=0;i<numberColumns;i++) {
1120        if (objective[i]) {
1121          column[n]=i;
1122          element[n++]=objective[i];
1123          objective[i]=0.0;
1124        }
1125      }
1126      column[n]=objectiveVariable_;
1127      element[n++]=-1.0;
1128      starts[1]=n;
1129      double offset = - modelPtr_->objectiveOffset();
1130      assert (!offset); // get sign right if happens
1131      modelPtr_->setObjectiveOffset(0.0);
1132      double lowerBound = -COIN_DBL_MAX;
1133      addRows(1,starts,column,element,&lowerBound,&offset);
1134      delete [] column;
1135      delete [] element;
1136    }
1137    OsiObject ** objects = new OsiObject * [nBi+nInt];
1138    char * marked = new char [numberColumns];
1139    memset(marked,0,numberColumns);
1140    // statistics I-I I-x x-x
1141    int stats[3]={0,0,0};
1142    double * sort = new double [nBi];
1143    nBi=nInt;
1144    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1145      if (quadraticObjective)
1146        startQuadratic[iColumn] = numberQuadratic;
1147      // See if quadratic objective
1148      const char * expr = coinModel_.getColumnObjectiveAsString(iColumn);
1149      if (strcmp(expr,"Numeric")) {
1150        // need bounds
1151        fakeBounds(this,iColumn,defaultBound_);
1152        // value*x*y
1153        char temp[20000];
1154        strcpy(temp,expr);
1155        char * pos = temp;
1156        bool ifFirst=true;
1157        while (*pos) {
1158          double value;
1159          int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1160          // must be column unless first when may be linear term
1161          if (jColumn>=0) {
1162            if (quadraticObjective) {
1163              columnQuadratic[numberQuadratic]=jColumn;
1164              elementQuadratic[numberQuadratic++]=2.0*value; // convention
1165            }
1166            // need bounds
1167            fakeBounds(this,jColumn,defaultBound_);
1168            double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
1169            if (meshI)
1170              marked[iColumn]=1;
1171            double meshJ = coinModel_.isInteger(jColumn) ? 1.0 : 0.0;
1172            if (meshJ)
1173              marked[jColumn]=1;
1174            // stats etc
1175            if (meshI) {
1176              if (meshJ)
1177                stats[0]++;
1178              else
1179                stats[1]++;
1180            } else {
1181              if (meshJ)
1182                stats[1]++;
1183              else
1184                stats[2]++;
1185            }
1186            if (iColumn<=jColumn)
1187              sort[nBi-nInt]=iColumn + numberColumns*jColumn;
1188            else
1189              sort[nBi-nInt]=jColumn + numberColumns*iColumn;
1190            if (!meshJ&&!meshI) {
1191              meshI=defaultMeshSize_;
1192              meshJ=0.0;
1193            }
1194            OsiBiLinear * newObj = new OsiBiLinear(this,iColumn,jColumn,objectiveRow_,value,meshI,meshJ,
1195                                                   nBi-nInt,(const OsiObject **) (objects+nInt));
1196            newObj->setPriority(biLinearPriority_);
1197            objects[nBi++] = newObj;
1198          } else if (jColumn==-2) {
1199          } else {
1200            printf("bad nonlinear term %s\n",temp);
1201            abort();
1202          }
1203          ifFirst=false;
1204        }
1205      }
1206    }
1207    // stats
1208    printf("There were %d I-I, %d I-x and %d x-x bilinear in objective\n",
1209           stats[0],stats[1],stats[2]);
1210    if (quadraticObjective) {
1211      startQuadratic[numberColumns] = numberQuadratic;
1212      quadraticModel_->loadQuadraticObjective(numberColumns,startQuadratic,
1213                                              columnQuadratic,elementQuadratic);
1214      delete [] startQuadratic;
1215      delete [] columnQuadratic;
1216      delete [] elementQuadratic;
1217    }
1218    for (iRow=0;iRow<numberRows;iRow++) {   
1219      CoinModelLink triple=coinModel_.firstInRow(iRow);
1220      while (triple.column()>=0) {
1221        int iColumn = triple.column();
1222        const char *  el = coinModel_.getElementAsString(iRow,iColumn);
1223        if (strcmp("Numeric",el)) {
1224          // need bounds
1225          fakeBounds(this,iColumn,defaultBound_);
1226          // value*x*y
1227          char temp[20000];
1228          strcpy(temp,el);
1229          char * pos = temp;
1230          bool ifFirst=true;
1231          while (*pos) {
1232            double value;
1233            int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
1234            // must be column unless first when may be linear term
1235            if (jColumn>=0) {
1236              // need bounds
1237              fakeBounds(this,jColumn,defaultBound_);
1238              double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
1239              if (meshI)
1240                marked[iColumn]=1;
1241              double meshJ = coinModel_.isInteger(jColumn) ? 1.0 : 0.0;
1242              if (meshJ)
1243                marked[jColumn]=1;
1244              // stats etc
1245              if (meshI) {
1246                if (meshJ)
1247                  stats[0]++;
1248                else
1249                  stats[1]++;
1250              } else {
1251                if (meshJ)
1252                  stats[1]++;
1253                else
1254                  stats[2]++;
1255              }
1256              if (iColumn<=jColumn)
1257                sort[nBi-nInt]=iColumn + numberColumns*jColumn;
1258              else
1259                sort[nBi-nInt]=jColumn + numberColumns*iColumn;
1260              if (!meshJ&&!meshI) {
1261                meshI=defaultMeshSize_;
1262                meshJ=0.0;
1263              }
1264              OsiBiLinear * newObj = new OsiBiLinear(this,iColumn,jColumn,iRow,value,meshI,meshJ,
1265                                                     nBi-nInt,(const OsiObject **) (objects+nInt));
1266              newObj->setPriority(biLinearPriority_);
1267              objects[nBi++] = newObj;
1268            } else if (jColumn==-2) {
1269            } else {
1270              printf("bad nonlinear term %s\n",temp);
1271              abort();
1272            }
1273            ifFirst=false;
1274          }
1275        }
1276        triple=coinModel_.next(triple);
1277      }
1278    }
1279    {
1280      // stats
1281      std::sort(sort,sort+nBi-nInt);
1282      int nDiff=0;
1283      double last=-1.0;
1284      for (int i=0;i<nBi-nInt;i++) {
1285        if (sort[i]!=last)
1286          nDiff++;
1287        last=sort[i];
1288      }
1289      delete [] sort;
1290      printf("There were %d I-I, %d I-x and %d x-x bilinear in total of which %d were duplicates\n",
1291             stats[0],stats[1],stats[2],nBi-nInt-nDiff);
1292    }
1293    nInt=0;
1294    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
1295      if (coinModel_.isInteger(iColumn)) {
1296        objects[nInt] = new OsiSimpleInteger(this,iColumn);
1297        if (marked[iColumn])
1298          objects[nInt]->setPriority(integerPriority_);
1299        else
1300          objects[nInt]->setPriority(integerPriority_);
1301        nInt++;
1302      }
1303    }
1304    nInt=nBi;
1305    delete [] marked;
1306    if (numberErrors) {
1307      // errors
1308      gutsOfDestructor();
1309      numberVariables_=-1;
1310    } else {
1311      addObjects(nInt,objects);
1312      int i;
1313      for (i=0;i<nInt;i++)
1314        delete objects[i];
1315      delete [] objects;
1316      // Now do dummy bound stuff
1317      matrix_ = new CoinPackedMatrix(*getMatrixByCol());
1318      info_ = new OsiLinkedBound [numberVariables_];
1319      for ( i=0;i<numberVariables_;i++) {
1320        info_[i] = OsiLinkedBound(this,which[i],0,NULL,NULL,NULL);
1321      }
1322      // Do row copy but just part
1323      int numberRows2 = objectiveRow_>=0 ? numberRows+1 : numberRows;
1324      int * whichRows = new int [numberRows2];
1325      int * whichColumns = new int [numberColumns];
1326      CoinIotaN(whichRows,numberRows2,0);
1327      CoinIotaN(whichColumns,numberColumns,0);
1328      originalRowCopy_ = new CoinPackedMatrix(*getMatrixByRow(),
1329                                              numberRows2,whichRows,
1330                                              numberColumns,whichColumns);
1331      delete [] whichColumns;
1332      numberNonLinearRows_=0;
1333      CoinZeroN(whichRows,numberRows2);
1334      for ( i =0;i<numberObjects_;i++) {
1335        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1336        if (obj) {
1337          int xyRow = obj->xyRow();
1338          assert (xyRow>=0&&xyRow<numberRows2); // even if obj we should move
1339          whichRows[xyRow]++;
1340        }
1341      }
1342      int * pos = new int [numberRows2];
1343      int n=0;
1344      for (i=0;i<numberRows2;i++) {
1345        if (whichRows[i]) {
1346          pos[numberNonLinearRows_]=n;
1347          n+=whichRows[i]; 
1348          whichRows[i]=numberNonLinearRows_;
1349          numberNonLinearRows_++;
1350        } else {
1351          whichRows[i]=-1;
1352        }
1353      }
1354      startNonLinear_ = new int [numberNonLinearRows_+1];
1355      memcpy(startNonLinear_,pos,numberNonLinearRows_*sizeof(int));
1356      startNonLinear_[numberNonLinearRows_]=n;
1357      rowNonLinear_ = new int [numberNonLinearRows_];
1358      convex_ = new int [numberNonLinearRows_];
1359      // do row numbers now
1360      numberNonLinearRows_=0;
1361      for (i=0;i<numberRows2;i++) {
1362        if (whichRows[i]>=0) {
1363          rowNonLinear_[numberNonLinearRows_++]=i;
1364        }
1365      }
1366      whichNonLinear_ = new int [n];
1367      for ( i =0;i<numberObjects_;i++) {
1368        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
1369        if (obj) {
1370          int xyRow = obj->xyRow();
1371          int k=whichRows[xyRow];
1372          int put = pos[k];
1373          pos[k]++;
1374          whichNonLinear_[put]=i;
1375        }
1376      }
1377      delete [] pos;
1378      delete [] whichRows;
1379      analyzeObjects();
1380    }
1381  }
1382  // See if there are any quadratic bounds
1383  int nQ=0;
1384  const CoinPackedMatrix * rowCopy = getMatrixByRow();
1385  //const double * element = rowCopy->getElements();
1386  //const int * column = rowCopy->getIndices();
1387  //const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1388  const int * rowLength = rowCopy->getVectorLengths();
1389  const double * rowLower = getRowLower();
1390  const double * rowUpper = getRowUpper();
1391  for (int iObject =0;iObject<numberObjects_;iObject++) {
1392    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[iObject]);
1393    if (obj) {
1394      int xyRow = obj->xyRow();
1395      if (rowLength[xyRow]==4) {
1396        // we have simple bound
1397        nQ++;
1398        double coefficient = obj->coefficient();
1399        double lo = rowLower[xyRow];
1400        double up = rowUpper[xyRow];
1401        if (coefficient!=1.0) {
1402          printf("*** double check code here\n");
1403          if (coefficient<0.0) {
1404            double temp = lo;
1405            lo = - up;
1406            up = - temp;
1407            coefficient = - coefficient;
1408          }
1409          if (lo>-1.0e20)
1410            lo /= coefficient;
1411          if (up <1.0e20)
1412            up /= coefficient;
1413          setRowLower(xyRow,lo);
1414          setRowUpper(xyRow,up);
1415          // we also need to change elements in matrix_
1416        }
1417        int type=0;
1418        if (lo==up) {
1419          // good news
1420          type=3;
1421          coefficient = lo;
1422        } else if (lo<-1.0e20) {
1423          assert (up<1.0e20);
1424          coefficient = up;
1425          type = 1;
1426          // can we make equality?
1427        } else if (up>1.0e20) {
1428          coefficient = lo;
1429          type = 2;
1430          // can we make equality?
1431        } else {
1432          // we would need extra code
1433          abort();
1434        }
1435        obj->setBoundType(type);
1436        obj->setCoefficient(coefficient);
1437        // can do better if integer?
1438        assert (!isInteger(obj->xColumn()));
1439        assert (!isInteger(obj->yColumn()));
1440      }
1441    }
1442  }
1443  delete [] which;
1444}
1445// Analyze constraints to see which are convex (quadratic)
1446void 
1447OsiSolverLink::analyzeObjects()
1448{
1449  // space for starts
1450  int numberColumns = coinModel_.numberColumns();
1451  int * start = new int [numberColumns+1];
1452  const double * rowLower = getRowLower();
1453  const double * rowUpper = getRowUpper();
1454  for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
1455    int iRow = rowNonLinear_[iNon];
1456    int numberElements = startNonLinear_[iNon+1]-startNonLinear_[iNon];
1457    // triplet arrays
1458    int * iColumn = new int [2*numberElements+1];
1459    int * jColumn = new int [2*numberElements];
1460    double * element = new double [2*numberElements];
1461    int i;
1462    int n=0;
1463    for ( i =startNonLinear_[iNon];i<startNonLinear_[iNon+1];i++) {
1464      OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
1465      assert (obj);
1466      int xColumn = obj->xColumn();
1467      int yColumn = obj->yColumn();
1468      double coefficient = obj->coefficient();
1469      if (xColumn!=yColumn) {
1470        iColumn[n]=xColumn;
1471        jColumn[n]=yColumn;
1472        element[n++]=coefficient;
1473        iColumn[n]=yColumn;
1474        jColumn[n]=xColumn;
1475        element[n++]=coefficient;
1476      } else {
1477        iColumn[n]=xColumn;
1478        jColumn[n]=xColumn;
1479        element[n++]=coefficient;
1480      }
1481    }
1482    // First sort in column order
1483    CoinSort_3(iColumn,iColumn+n,jColumn,element);
1484    // marker at end
1485    iColumn[n]=numberColumns;
1486    int lastI=iColumn[0];
1487    // compute starts
1488    start[0]=0;
1489    for (i=1;i<n+1;i++) {
1490      if (iColumn[i]!=lastI) {
1491        while (lastI<iColumn[i]) {
1492          start[lastI+1]=i;
1493          lastI++;
1494        }
1495        lastI=iColumn[i];
1496      }
1497    }
1498    // -1 unknown, 0 convex, 1 nonconvex
1499    int status=-1;
1500    int statusNegative=-1;
1501    int numberLong=0; // number with >2 elements
1502    for (int k=0;k<numberColumns;k++) {
1503      int first = start[k];
1504      int last = start[k+1];
1505      if (last>first) {
1506        int j;
1507        double diagonal = 0.0;
1508        int whichK=-1;
1509        for (j=first;j<last;j++) {
1510          if (jColumn[j]==k) {
1511            diagonal = element[j];
1512            status = diagonal >0 ? 0 : 1;
1513            statusNegative = diagonal <0 ? 0 : 1;
1514            whichK = (j==first) ? j+1 :j-1;
1515            break;
1516          }
1517        }
1518        if (last==first+1) {
1519          // just one entry
1520          if (!diagonal) {
1521            // one off diagonal - not positive semi definite
1522            status=1;
1523            statusNegative=1;
1524          }
1525        } else if (diagonal) {
1526          if (last==first+2) {
1527            // other column and element
1528            double otherElement=element[whichK];;
1529            int otherColumn = jColumn[whichK];
1530            double otherDiagonal=0.0;
1531            // check 2x2 determinant - unless past and 2 long
1532            if (otherColumn>i||start[otherColumn+1]>start[otherColumn]+2) {
1533              for (j=start[otherColumn];j<start[otherColumn+1];j++) {
1534                if (jColumn[j]==otherColumn) {
1535                  otherDiagonal = element[j];
1536                  break;
1537                }
1538              }
1539              // determinant
1540              double determinant = diagonal*otherDiagonal-otherElement*otherElement;
1541              if (determinant<-1.0e-12) {
1542                // not positive semi definite
1543                status=1;
1544                statusNegative=1;
1545              } else if (start[otherColumn+1]>start[otherColumn]+2&&determinant<1.0e-12) {
1546                // not positive semi definite
1547                status=1;
1548                statusNegative=1;
1549              }
1550            }
1551          } else {
1552            numberLong++;
1553          }
1554        }
1555      }
1556    }
1557    if ((status==0||statusNegative==0)&&numberLong) {
1558      // need to do more work
1559      printf("Needs more work\n");
1560    }
1561    assert (status>0||statusNegative>0);
1562    if (!status) {
1563      convex_[iNon]=1;
1564      // equality may be ok
1565      assert (rowUpper[iRow]<1.0e20);
1566      specialOptions2_ |= 8;
1567    } else if (!statusNegative) {
1568      convex_[iNon]=-1;
1569      // equality may be ok
1570      assert (rowLower[iRow]>-1.0e20);
1571      specialOptions2_ |= 8;
1572    } else {
1573      convex_[iNon]=0;
1574    }
1575    printf("Convexity of row %d is %d\n",iRow,convex_[iNon]);
1576    delete [] iColumn;
1577    delete [] jColumn;
1578    delete [] element;
1579  }
1580  delete [] start;
1581}
1582//-------------------------------------------------------------------
1583// Clone
1584//-------------------------------------------------------------------
1585OsiSolverInterface * 
1586OsiSolverLink::clone(bool copyData) const
1587{
1588  assert (copyData);
1589  return new OsiSolverLink(*this);
1590}
1591
1592
1593//-------------------------------------------------------------------
1594// Copy constructor
1595//-------------------------------------------------------------------
1596OsiSolverLink::OsiSolverLink (
1597                  const OsiSolverLink & rhs)
1598  : OsiClpSolverInterface(rhs)
1599{
1600  gutsOfDestructor(true);
1601  gutsOfCopy(rhs);
1602  // something odd happens - try this
1603  OsiSolverInterface::operator=(rhs);
1604}
1605
1606//-------------------------------------------------------------------
1607// Destructor
1608//-------------------------------------------------------------------
1609OsiSolverLink::~OsiSolverLink ()
1610{
1611  gutsOfDestructor();
1612}
1613
1614//-------------------------------------------------------------------
1615// Assignment operator
1616//-------------------------------------------------------------------
1617OsiSolverLink &
1618OsiSolverLink::operator=(const OsiSolverLink& rhs)
1619{
1620  if (this != &rhs) { 
1621    gutsOfDestructor();
1622    OsiClpSolverInterface::operator=(rhs);
1623    gutsOfCopy(rhs);
1624  }
1625  return *this;
1626}
1627void 
1628OsiSolverLink::gutsOfDestructor(bool justNullify)
1629{
1630  if (!justNullify) {
1631    delete matrix_;
1632    delete originalRowCopy_;
1633    delete [] info_;
1634    delete [] bestSolution_;
1635    delete quadraticModel_;
1636    delete [] startNonLinear_;
1637    delete [] rowNonLinear_;
1638    delete [] convex_;
1639    delete [] whichNonLinear_;
1640  } 
1641  matrix_ = NULL;
1642  originalRowCopy_ = NULL;
1643  quadraticModel_ = NULL;
1644  numberNonLinearRows_=0;
1645  startNonLinear_ = NULL;
1646  rowNonLinear_ = NULL;
1647  convex_ = NULL;
1648  whichNonLinear_ = NULL;
1649  cbcModel_ = NULL;
1650  info_ = NULL;
1651  numberVariables_ = 0;
1652  specialOptions2_ = 0;
1653  objectiveRow_=-1;
1654  objectiveVariable_=-1;
1655  bestSolution_ = NULL;
1656  bestObjectiveValue_ =1.0e100;
1657  defaultMeshSize_ = 1.0e-4;
1658  defaultBound_ = 1.0e4;
1659  integerPriority_ = 1000;
1660  biLinearPriority_ = 10000;
1661}
1662void 
1663OsiSolverLink::gutsOfCopy(const OsiSolverLink & rhs)
1664{
1665  coinModel_ = rhs.coinModel_;
1666  numberVariables_ = rhs.numberVariables_;
1667  numberNonLinearRows_ = rhs.numberNonLinearRows_;
1668  specialOptions2_ = rhs.specialOptions2_;
1669  objectiveRow_=rhs.objectiveRow_;
1670  objectiveVariable_=rhs.objectiveVariable_;
1671  bestObjectiveValue_ = rhs.bestObjectiveValue_;
1672  defaultMeshSize_ = rhs.defaultMeshSize_;
1673  defaultBound_ = rhs.defaultBound_;
1674  integerPriority_ = rhs.integerPriority_;
1675  biLinearPriority_ = rhs.biLinearPriority_;
1676  cbcModel_ = rhs.cbcModel_;
1677  if (numberVariables_) { 
1678    if (rhs.matrix_)
1679      matrix_ = new CoinPackedMatrix(*rhs.matrix_);
1680    else
1681      matrix_=NULL;
1682    if (rhs.originalRowCopy_)
1683      originalRowCopy_ = new CoinPackedMatrix(*rhs.originalRowCopy_);
1684    else
1685      originalRowCopy_=NULL;
1686    info_ = new OsiLinkedBound [numberVariables_];
1687    for (int i=0;i<numberVariables_;i++) {
1688      info_[i] = OsiLinkedBound(rhs.info_[i]);
1689    }
1690    if (rhs.bestSolution_) {
1691      bestSolution_ = CoinCopyOfArray(rhs.bestSolution_,modelPtr_->getNumCols());
1692    } else {
1693      bestSolution_=NULL;
1694    }
1695  }
1696  if (numberNonLinearRows_) {
1697    startNonLinear_ = CoinCopyOfArray(rhs.startNonLinear_,numberNonLinearRows_+1);
1698    rowNonLinear_ = CoinCopyOfArray(rhs.rowNonLinear_,numberNonLinearRows_);
1699    convex_ = CoinCopyOfArray(rhs.convex_,numberNonLinearRows_);
1700    int numberEntries = startNonLinear_[numberNonLinearRows_];
1701    whichNonLinear_ = CoinCopyOfArray(rhs.whichNonLinear_,numberEntries);
1702  }
1703  if (rhs.quadraticModel_) {
1704    quadraticModel_ = new ClpSimplex(*rhs.quadraticModel_);
1705  } else {
1706    quadraticModel_ = NULL;
1707  }
1708}
1709// Add a bound modifier
1710void 
1711OsiSolverLink::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, int whichVariableAffected, 
1712                                double multiplier)
1713{
1714  bool found=false;
1715  int i;
1716  for ( i=0;i<numberVariables_;i++) {
1717    if (info_[i].variable()==whichVariable) {
1718      found=true;
1719      break;
1720    }
1721  }
1722  if (!found) {
1723    // add in
1724    OsiLinkedBound * temp = new OsiLinkedBound [numberVariables_+1];
1725    for (int i=0;i<numberVariables_;i++) 
1726      temp[i]= info_[i];
1727    delete [] info_;
1728    info_=temp;
1729    info_[numberVariables_++] = OsiLinkedBound(this,whichVariable,0,NULL,NULL,NULL);
1730  }
1731  info_[i].addBoundModifier(upperBoundAffected, useUpperBound,whichVariableAffected,multiplier);
1732}
1733// Update coefficients
1734int
1735OsiSolverLink::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
1736{
1737  double * lower = solver->columnLower();
1738  double * upper = solver->columnUpper();
1739  double * objective = solver->objective();
1740  int numberChanged=0;
1741  for (int iObject =0;iObject<numberObjects_;iObject++) {
1742    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[iObject]);
1743    if (obj) {
1744      numberChanged +=obj->updateCoefficients(lower,upper,objective,matrix,&basis_);
1745    }
1746  }
1747  return numberChanged;
1748}
1749// Set best solution found internally
1750void 
1751OsiSolverLink::setBestSolution(const double * solution, int numberColumns)
1752{
1753  delete [] bestSolution_;
1754  int numberColumnsThis = modelPtr_->numberColumns();
1755  bestSolution_ = new double [numberColumnsThis];
1756  CoinZeroN(bestSolution_,numberColumnsThis);
1757  memcpy(bestSolution_,solution,CoinMin(numberColumns,numberColumnsThis)*sizeof(double));
1758}
1759//#############################################################################
1760// Constructors, destructors  and assignment
1761//#############################################################################
1762
1763//-------------------------------------------------------------------
1764// Default Constructor
1765//-------------------------------------------------------------------
1766OsiLinkedBound::OsiLinkedBound ()
1767{
1768  model_ = NULL;
1769  variable_ = -1;
1770  numberAffected_ = 0;
1771  maximumAffected_ = numberAffected_;
1772  affected_ = NULL;
1773}
1774// Useful Constructor
1775OsiLinkedBound::OsiLinkedBound(OsiSolverInterface * model, int variable,
1776                               int numberAffected, const int * positionL, 
1777                               const int * positionU, const double * multiplier)
1778{
1779  model_ = model;
1780  variable_ = variable;
1781  numberAffected_ = 2*numberAffected;
1782  maximumAffected_ = numberAffected_;
1783  if (numberAffected_) { 
1784    affected_ = new boundElementAction[numberAffected_];
1785    int n=0;
1786    for (int i=0;i<numberAffected;i++) {
1787      // LB
1788      boundElementAction action;
1789      action.affect=2;
1790      action.ubUsed=0;
1791      action.type=0;
1792      action.affected=positionL[i];
1793      action.multiplier=multiplier[i];
1794      affected_[n++]=action;
1795      // UB
1796      action.affect=2;
1797      action.ubUsed=1;
1798      action.type=0;
1799      action.affected=positionU[i];
1800      action.multiplier=multiplier[i];
1801      affected_[n++]=action;
1802    }
1803  } else {
1804    affected_ = NULL;
1805  }
1806}
1807
1808//-------------------------------------------------------------------
1809// Copy constructor
1810//-------------------------------------------------------------------
1811OsiLinkedBound::OsiLinkedBound (
1812                  const OsiLinkedBound & rhs)
1813{
1814  model_ = rhs.model_;
1815  variable_ = rhs.variable_;
1816  numberAffected_ = rhs.numberAffected_;
1817  maximumAffected_ = rhs.maximumAffected_;
1818  if (numberAffected_) { 
1819    affected_ = new boundElementAction[maximumAffected_];
1820    memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
1821  } else {
1822    affected_ = NULL;
1823  }
1824}
1825
1826//-------------------------------------------------------------------
1827// Destructor
1828//-------------------------------------------------------------------
1829OsiLinkedBound::~OsiLinkedBound ()
1830{
1831  delete [] affected_;
1832}
1833
1834//-------------------------------------------------------------------
1835// Assignment operator
1836//-------------------------------------------------------------------
1837OsiLinkedBound &
1838OsiLinkedBound::operator=(const OsiLinkedBound& rhs)
1839{
1840  if (this != &rhs) { 
1841    delete [] affected_;
1842    model_ = rhs.model_;
1843    variable_ = rhs.variable_;
1844    numberAffected_ = rhs.numberAffected_;
1845    maximumAffected_ = rhs.maximumAffected_;
1846    if (numberAffected_) { 
1847      affected_ = new boundElementAction[maximumAffected_];
1848      memcpy(affected_,rhs.affected_,numberAffected_*sizeof(boundElementAction));
1849    } else {
1850      affected_ = NULL;
1851    }
1852  }
1853  return *this;
1854}
1855// Add a bound modifier
1856void 
1857OsiLinkedBound::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, 
1858                                 double multiplier)
1859{
1860  if (numberAffected_==maximumAffected_) {
1861    maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
1862    boundElementAction * temp = new boundElementAction[maximumAffected_];
1863    memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
1864    delete [] affected_;
1865    affected_ = temp;
1866  }
1867  boundElementAction action;
1868  action.affect=upperBoundAffected ? 1 : 0;
1869  action.ubUsed=useUpperBound ? 1 : 0;
1870  action.type=2;
1871  action.affected=whichVariable;
1872  action.multiplier=multiplier;
1873  affected_[numberAffected_++]=action;
1874 
1875}
1876// Update other bounds
1877void 
1878OsiLinkedBound::updateBounds(ClpSimplex * solver)
1879{
1880  double * lower = solver->columnLower();
1881  double * upper = solver->columnUpper();
1882  double lo = lower[variable_];
1883  double up = upper[variable_];
1884  // printf("bounds for %d are %g and %g\n",variable_,lo,up);
1885  for (int j=0;j<numberAffected_;j++) {
1886    if (affected_[j].affect<2) {
1887      double multiplier = affected_[j].multiplier;
1888      assert (affected_[j].type==2);
1889      int iColumn = affected_[j].affected;
1890      double useValue = (affected_[j].ubUsed) ? up : lo;
1891      if (affected_[j].affect==0) 
1892        lower[iColumn] = CoinMin(upper[iColumn],CoinMax(lower[iColumn],multiplier*useValue));
1893      else
1894        upper[iColumn] = CoinMax(lower[iColumn],CoinMin(upper[iColumn],multiplier*useValue));
1895    }
1896  }
1897}
1898#if 0
1899// Add an element modifier
1900void
1901OsiLinkedBound::addCoefficientModifier(bool useUpperBound, int position,
1902                                 double multiplier)
1903{
1904  if (numberAffected_==maximumAffected_) {
1905    maximumAffected_ = maximumAffected_+10+maximumAffected_/4;
1906    boundElementAction * temp = new boundElementAction[maximumAffected_];
1907    memcpy(temp,affected_,numberAffected_*sizeof(boundElementAction));
1908    delete [] affected_;
1909    affected_ = temp;
1910  }
1911  boundElementAction action;
1912  action.affect=2;
1913  action.ubUsed=useUpperBound ? 1 : 0;
1914  action.type=0;
1915  action.affected=position;
1916  action.multiplier=multiplier;
1917  affected_[numberAffected_++]=action;
1918 
1919}
1920// Update coefficients
1921void
1922OsiLinkedBound::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
1923{
1924  double * lower = solver->columnLower();
1925  double * upper = solver->columnUpper();
1926  double * element = matrix->getMutableElements();
1927  double lo = lower[variable_];
1928  double up = upper[variable_];
1929  // printf("bounds for %d are %g and %g\n",variable_,lo,up);
1930  for (int j=0;j<numberAffected_;j++) {
1931    if (affected_[j].affect==2) {
1932      double multiplier = affected_[j].multiplier;
1933      assert (affected_[j].type==0);
1934      int position = affected_[j].affected;
1935      //double old = element[position];
1936      if (affected_[j].ubUsed)
1937        element[position] = multiplier*up;
1938      else
1939        element[position] = multiplier*lo;
1940      //if ( old != element[position])
1941      //printf("change at %d from %g to %g\n",position,old,element[position]);
1942    }
1943  }
1944}
1945#endif
1946// Default Constructor
1947CbcHeuristicDynamic3::CbcHeuristicDynamic3() 
1948  :CbcHeuristic()
1949{
1950}
1951
1952// Constructor from model
1953CbcHeuristicDynamic3::CbcHeuristicDynamic3(CbcModel & model)
1954  :CbcHeuristic(model)
1955{
1956}
1957
1958// Destructor
1959CbcHeuristicDynamic3::~CbcHeuristicDynamic3 ()
1960{
1961}
1962
1963// Clone
1964CbcHeuristic *
1965CbcHeuristicDynamic3::clone() const
1966{
1967  return new CbcHeuristicDynamic3(*this);
1968}
1969
1970// Copy constructor
1971CbcHeuristicDynamic3::CbcHeuristicDynamic3(const CbcHeuristicDynamic3 & rhs)
1972:
1973  CbcHeuristic(rhs)
1974{
1975}
1976
1977// Returns 1 if solution, 0 if not
1978int
1979CbcHeuristicDynamic3::solution(double & solutionValue,
1980                         double * betterSolution)
1981{
1982  if (!model_)
1983    return 0;
1984  OsiSolverLink * clpSolver
1985    = dynamic_cast<OsiSolverLink *> (model_->solver());
1986  assert (clpSolver); 
1987  double newSolutionValue = clpSolver->bestObjectiveValue();
1988  const double * solution = clpSolver->bestSolution();
1989  if (newSolutionValue<solutionValue&&solution) {
1990    int numberColumns = clpSolver->getNumCols();
1991    // new solution
1992    memcpy(betterSolution,solution,numberColumns*sizeof(double));
1993    solutionValue = newSolutionValue;
1994    return 1;
1995  } else {
1996    return 0;
1997  }
1998}
1999// update model
2000void CbcHeuristicDynamic3::setModel(CbcModel * model)
2001{
2002  model_ = model;
2003}
2004// Resets stuff if model changes
2005void 
2006CbcHeuristicDynamic3::resetModel(CbcModel * model)
2007{
2008  model_ = model;
2009}
2010#include <cassert>
2011#include <cmath>
2012#include <cfloat>
2013//#define CBC_DEBUG
2014
2015#include "OsiSolverInterface.hpp"
2016//#include "OsiBranchLink.hpp"
2017#include "CoinError.hpp"
2018#include "CoinHelperFunctions.hpp"
2019#include "CoinPackedMatrix.hpp"
2020#include "CoinWarmStartBasis.hpp"
2021
2022// Default Constructor
2023OsiOldLink::OsiOldLink ()
2024  : OsiSOS(),
2025    numberLinks_(0)
2026{
2027}
2028
2029// Useful constructor (which are indices)
2030OsiOldLink::OsiOldLink (const OsiSolverInterface * solver,  int numberMembers,
2031           int numberLinks, int first , const double * weights, int identifier)
2032  : OsiSOS(),
2033    numberLinks_(numberLinks)
2034{
2035  numberMembers_ = numberMembers;
2036  members_ = NULL;
2037  sosType_ = 1;
2038  if (numberMembers_) {
2039    weights_ = new double[numberMembers_];
2040    members_ = new int[numberMembers_*numberLinks_];
2041    if (weights) {
2042      memcpy(weights_,weights,numberMembers_*sizeof(double));
2043    } else {
2044      for (int i=0;i<numberMembers_;i++)
2045        weights_[i]=i;
2046    }
2047    // weights must be increasing
2048    int i;
2049    double last=-COIN_DBL_MAX;
2050    for (i=0;i<numberMembers_;i++) {
2051      assert (weights_[i]>last+1.0e-12);
2052      last=weights_[i];
2053    }
2054    for (i=0;i<numberMembers_*numberLinks_;i++) {
2055      members_[i]=first+i;
2056    }
2057  } else {
2058    weights_ = NULL;
2059  }
2060}
2061
2062// Useful constructor (which are indices)
2063OsiOldLink::OsiOldLink (const OsiSolverInterface * solver,  int numberMembers,
2064           int numberLinks, int sosType, const int * which , const double * weights, int identifier)
2065  : OsiSOS(),
2066    numberLinks_(numberLinks)
2067{
2068  numberMembers_ = numberMembers;
2069  members_ = NULL;
2070  sosType_ = 1;
2071  if (numberMembers_) {
2072    weights_ = new double[numberMembers_];
2073    members_ = new int[numberMembers_*numberLinks_];
2074    if (weights) {
2075      memcpy(weights_,weights,numberMembers_*sizeof(double));
2076    } else {
2077      for (int i=0;i<numberMembers_;i++)
2078        weights_[i]=i;
2079    }
2080    // weights must be increasing
2081    int i;
2082    double last=-COIN_DBL_MAX;
2083    for (i=0;i<numberMembers_;i++) {
2084      assert (weights_[i]>last+1.0e-12);
2085      last=weights_[i];
2086    }
2087    for (i=0;i<numberMembers_*numberLinks_;i++) {
2088      members_[i]= which[i];
2089    }
2090  } else {
2091    weights_ = NULL;
2092  }
2093}
2094
2095// Copy constructor
2096OsiOldLink::OsiOldLink ( const OsiOldLink & rhs)
2097  :OsiSOS(rhs)
2098{
2099  numberLinks_ = rhs.numberLinks_;
2100  if (numberMembers_) {
2101    delete [] members_;
2102    members_ = CoinCopyOfArray(rhs.members_,numberMembers_*numberLinks_);
2103  }
2104}
2105
2106// Clone
2107OsiObject *
2108OsiOldLink::clone() const
2109{
2110  return new OsiOldLink(*this);
2111}
2112
2113// Assignment operator
2114OsiOldLink & 
2115OsiOldLink::operator=( const OsiOldLink& rhs)
2116{
2117  if (this!=&rhs) {
2118    OsiSOS::operator=(rhs);
2119    delete [] members_;
2120    numberLinks_ = rhs.numberLinks_;
2121    if (numberMembers_) {
2122      members_ = CoinCopyOfArray(rhs.members_,numberMembers_*numberLinks_);
2123    } else {
2124      members_ = NULL;
2125    }
2126  }
2127  return *this;
2128}
2129
2130// Destructor
2131OsiOldLink::~OsiOldLink ()
2132{
2133}
2134
2135// Infeasibility - large is 0.5
2136double 
2137OsiOldLink::infeasibility(const OsiBranchingInformation * info,int & whichWay) const
2138{
2139  int j;
2140  int firstNonZero=-1;
2141  int lastNonZero = -1;
2142  const double * solution = info->solution_;
2143  //const double * lower = info->lower_;
2144  const double * upper = info->upper_;
2145  double integerTolerance = info->integerTolerance_;
2146  double weight = 0.0;
2147  double sum =0.0;
2148
2149  // check bounds etc
2150  double lastWeight=-1.0e100;
2151  int base=0;
2152  for (j=0;j<numberMembers_;j++) {
2153    for (int k=0;k<numberLinks_;k++) {
2154      int iColumn = members_[base+k];
2155      if (lastWeight>=weights_[j]-1.0e-7)
2156        throw CoinError("Weights too close together in OsiLink","infeasibility","OsiLink");
2157      lastWeight = weights_[j];
2158      double value = CoinMax(0.0,solution[iColumn]);
2159      sum += value;
2160      if (value>integerTolerance&&upper[iColumn]) {
2161        // Possibly due to scaling a fixed variable might slip through
2162        if (value>upper[iColumn]+1.0e-8) {
2163#ifdef OSI_DEBUG
2164          printf("** Variable %d (%d) has value %g and upper bound of %g\n",
2165                 iColumn,j,value,upper[iColumn]);
2166#endif
2167        } 
2168        value = CoinMin(value,upper[iColumn]);
2169        weight += weights_[j]*value;
2170        if (firstNonZero<0)
2171          firstNonZero=j;
2172        lastNonZero=j;
2173      }
2174    }
2175    base += numberLinks_;
2176  }
2177  double valueInfeasibility;
2178  whichWay=1;
2179  whichWay_=1;
2180  if (lastNonZero-firstNonZero>=sosType_) {
2181    // find where to branch
2182    assert (sum>0.0);
2183    weight /= sum;
2184    valueInfeasibility = lastNonZero-firstNonZero+1;
2185    valueInfeasibility *= 0.5/((double) numberMembers_);
2186    //#define DISTANCE
2187#ifdef DISTANCE
2188    assert (sosType_==1); // code up
2189    /* may still be satisfied.
2190       For LOS type 2 we might wish to move coding around
2191       and keep initial info in model_ for speed
2192    */
2193    int iWhere;
2194    bool possible=false;
2195    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
2196      if (fabs(weight-weights_[iWhere])<1.0e-8) {
2197        possible=true;
2198        break;
2199      }
2200    }
2201    if (possible) {
2202      // One could move some of this (+ arrays) into model_
2203      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
2204      const double * element = matrix->getMutableElements();
2205      const int * row = matrix->getIndices();
2206      const CoinBigIndex * columnStart = matrix->getVectorStarts();
2207      const int * columnLength = matrix->getVectorLengths();
2208      const double * rowSolution = solver->getRowActivity();
2209      const double * rowLower = solver->getRowLower();
2210      const double * rowUpper = solver->getRowUpper();
2211      int numberRows = matrix->getNumRows();
2212      double * array = new double [numberRows];
2213      CoinZeroN(array,numberRows);
2214      int * which = new int [numberRows];
2215      int n=0;
2216      int base=numberLinks_*firstNonZero;
2217      for (j=firstNonZero;j<=lastNonZero;j++) {
2218        for (int k=0;k<numberLinks_;k++) {
2219          int iColumn = members_[base+k];
2220          double value = CoinMax(0.0,solution[iColumn]);
2221          if (value>integerTolerance&&upper[iColumn]) {
2222            value = CoinMin(value,upper[iColumn]);
2223            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
2224              int iRow = row[j];
2225              double a = array[iRow];
2226              if (a) {
2227                a += value*element[j];
2228                if (!a)
2229                  a = 1.0e-100;
2230              } else {
2231                which[n++]=iRow;
2232                a=value*element[j];
2233                assert (a);
2234              }
2235              array[iRow]=a;
2236            }
2237          }
2238        }
2239        base += numberLinks_;
2240      }
2241      base=numberLinks_*iWhere;
2242      for (int k=0;k<numberLinks_;k++) {
2243        int iColumn = members_[base+k];
2244        const double value = 1.0;
2245        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
2246          int iRow = row[j];
2247          double a = array[iRow];
2248          if (a) {
2249            a -= value*element[j];
2250            if (!a)
2251              a = 1.0e-100;
2252          } else {
2253            which[n++]=iRow;
2254            a=-value*element[j];
2255            assert (a);
2256          }
2257          array[iRow]=a;
2258        }
2259      }
2260      for (j=0;j<n;j++) {
2261        int iRow = which[j];
2262        // moving to point will increase row solution by this
2263        double distance = array[iRow];
2264        if (distance>1.0e-8) {
2265          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
2266            possible=false;
2267            break;
2268          }
2269        } else if (distance<-1.0e-8) {
2270          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
2271            possible=false;
2272            break;
2273          } 
2274        }
2275      }
2276      for (j=0;j<n;j++)
2277        array[which[j]]=0.0;
2278      delete [] array;
2279      delete [] which;
2280      if (possible) {
2281        valueInfeasibility=0.0;
2282        printf("possible %d %d %d\n",firstNonZero,lastNonZero,iWhere);
2283      }
2284    }
2285#endif
2286  } else {
2287    valueInfeasibility = 0.0; // satisfied
2288  }
2289  infeasibility_=valueInfeasibility;
2290  otherInfeasibility_=1.0-valueInfeasibility;
2291  return valueInfeasibility;
2292}
2293
2294// This looks at solution and sets bounds to contain solution
2295double
2296OsiOldLink::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
2297{
2298  int j;
2299  int firstNonZero=-1;
2300  int lastNonZero = -1;
2301  const double * solution = info->solution_;
2302  const double * upper = info->upper_;
2303  double integerTolerance = info->integerTolerance_;
2304  double weight = 0.0;
2305  double sum =0.0;
2306
2307  int base=0;
2308  for (j=0;j<numberMembers_;j++) {
2309    for (int k=0;k<numberLinks_;k++) {
2310      int iColumn = members_[base+k];
2311      double value = CoinMax(0.0,solution[iColumn]);
2312      sum += value;
2313      if (value>integerTolerance&&upper[iColumn]) {
2314        weight += weights_[j]*value;
2315        if (firstNonZero<0)
2316          firstNonZero=j;
2317        lastNonZero=j;
2318      }
2319    }
2320    base += numberLinks_;
2321  }
2322#ifdef DISTANCE
2323  if (lastNonZero-firstNonZero>sosType_-1) {
2324    /* may still be satisfied.
2325       For LOS type 2 we might wish to move coding around
2326       and keep initial info in model_ for speed
2327    */
2328    int iWhere;
2329    bool possible=false;
2330    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
2331      if (fabs(weight-weights_[iWhere])<1.0e-8) {
2332        possible=true;
2333        break;
2334      }
2335    }
2336    if (possible) {
2337      // One could move some of this (+ arrays) into model_
2338      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
2339      const double * element = matrix->getMutableElements();
2340      const int * row = matrix->getIndices();
2341      const CoinBigIndex * columnStart = matrix->getVectorStarts();
2342      const int * columnLength = matrix->getVectorLengths();
2343      const double * rowSolution = solver->getRowActivity();
2344      const double * rowLower = solver->getRowLower();
2345      const double * rowUpper = solver->getRowUpper();
2346      int numberRows = matrix->getNumRows();
2347      double * array = new double [numberRows];
2348      CoinZeroN(array,numberRows);
2349      int * which = new int [numberRows];
2350      int n=0;
2351      int base=numberLinks_*firstNonZero;
2352      for (j=firstNonZero;j<=lastNonZero;j++) {
2353        for (int k=0;k<numberLinks_;k++) {
2354          int iColumn = members_[base+k];
2355          double value = CoinMax(0.0,solution[iColumn]);
2356          if (value>integerTolerance&&upper[iColumn]) {
2357            value = CoinMin(value,upper[iColumn]);
2358            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
2359              int iRow = row[j];
2360              double a = array[iRow];
2361              if (a) {
2362                a += value*element[j];
2363                if (!a)
2364                  a = 1.0e-100;
2365              } else {
2366                which[n++]=iRow;
2367                a=value*element[j];
2368                assert (a);
2369              }
2370              array[iRow]=a;
2371            }
2372          }
2373        }
2374        base += numberLinks_;
2375      }
2376      base=numberLinks_*iWhere;
2377      for (int k=0;k<numberLinks_;k++) {
2378        int iColumn = members_[base+k];
2379        const double value = 1.0;
2380        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
2381          int iRow = row[j];
2382          double a = array[iRow];
2383          if (a) {
2384            a -= value*element[j];
2385            if (!a)
2386              a = 1.0e-100;
2387          } else {
2388            which[n++]=iRow;
2389            a=-value*element[j];
2390            assert (a);
2391          }
2392          array[iRow]=a;
2393        }
2394      }
2395      for (j=0;j<n;j++) {
2396        int iRow = which[j];
2397        // moving to point will increase row solution by this
2398        double distance = array[iRow];
2399        if (distance>1.0e-8) {
2400          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
2401            possible=false;
2402            break;
2403          }
2404        } else if (distance<-1.0e-8) {
2405          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
2406            possible=false;
2407            break;
2408          } 
2409        }
2410      }
2411      for (j=0;j<n;j++)
2412        array[which[j]]=0.0;
2413      delete [] array;
2414      delete [] which;
2415      if (possible) {
2416        printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere);
2417        firstNonZero=iWhere;
2418        lastNonZero=iWhere;
2419      }
2420    }
2421  }
2422#else
2423  assert (lastNonZero-firstNonZero<sosType_) ;
2424#endif
2425  base=0;
2426  for (j=0;j<firstNonZero;j++) {
2427    for (int k=0;k<numberLinks_;k++) {
2428      int iColumn = members_[base+k];
2429      solver->setColUpper(iColumn,0.0);
2430    }
2431    base += numberLinks_;
2432  }
2433  // skip
2434  base += numberLinks_;
2435  for (j=lastNonZero+1;j<numberMembers_;j++) {
2436    for (int k=0;k<numberLinks_;k++) {
2437      int iColumn = members_[base+k];
2438      solver->setColUpper(iColumn,0.0);
2439    }
2440    base += numberLinks_;
2441  }
2442  // go to coding as in OsiSOS
2443  abort();
2444  return -1.0;
2445}
2446
2447// Redoes data when sequence numbers change
2448void 
2449OsiOldLink::resetSequenceEtc(int numberColumns, const int * originalColumns)
2450{
2451  int n2=0;
2452  for (int j=0;j<numberMembers_*numberLinks_;j++) {
2453    int iColumn = members_[j];
2454    int i;
2455    for (i=0;i<numberColumns;i++) {
2456      if (originalColumns[i]==iColumn)
2457        break;
2458    }
2459    if (i<numberColumns) {
2460      members_[n2]=i;
2461      weights_[n2++]=weights_[j];
2462    }
2463  }
2464  if (n2<numberMembers_) {
2465    printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2/numberLinks_);
2466    numberMembers_=n2/numberLinks_;
2467  }
2468}
2469
2470// Creates a branching object
2471OsiBranchingObject * 
2472OsiOldLink::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
2473{
2474  int j;
2475  const double * solution = info->solution_;
2476  double tolerance = info->primalTolerance_;
2477  const double * upper = info->upper_;
2478  int firstNonFixed=-1;
2479  int lastNonFixed=-1;
2480  int firstNonZero=-1;
2481  int lastNonZero = -1;
2482  double weight = 0.0;
2483  double sum =0.0;
2484  int base=0;
2485  for (j=0;j<numberMembers_;j++) {
2486    for (int k=0;k<numberLinks_;k++) {
2487      int iColumn = members_[base+k];
2488      if (upper[iColumn]) {
2489        double value = CoinMax(0.0,solution[iColumn]);
2490        sum += value;
2491        if (firstNonFixed<0)
2492          firstNonFixed=j;
2493        lastNonFixed=j;
2494        if (value>tolerance) {
2495          weight += weights_[j]*value;
2496          if (firstNonZero<0)
2497            firstNonZero=j;
2498          lastNonZero=j;
2499        }
2500      }
2501    }
2502    base += numberLinks_;
2503  }
2504  assert (lastNonZero-firstNonZero>=sosType_) ;
2505  // find where to branch
2506  assert (sum>0.0);
2507  weight /= sum;
2508  int iWhere;
2509  double separator=0.0;
2510  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
2511    if (weight<weights_[iWhere+1])
2512      break;
2513  if (sosType_==1) {
2514    // SOS 1
2515    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
2516  } else {
2517    // SOS 2
2518    if (iWhere==firstNonFixed)
2519      iWhere++;;
2520    if (iWhere==lastNonFixed-1)
2521      iWhere = lastNonFixed-2;
2522    separator = weights_[iWhere+1];
2523  }
2524  // create object
2525  OsiBranchingObject * branch;
2526  branch = new OsiOldLinkBranchingObject(solver,this,way,separator);
2527  return branch;
2528}
2529OsiOldLinkBranchingObject::OsiOldLinkBranchingObject()
2530  :OsiSOSBranchingObject()
2531{
2532}
2533
2534// Useful constructor
2535OsiOldLinkBranchingObject::OsiOldLinkBranchingObject (OsiSolverInterface * solver,
2536                                              const OsiOldLink * set,
2537                                              int way ,
2538                                              double separator)
2539  :OsiSOSBranchingObject(solver,set,way,separator)
2540{
2541}
2542
2543// Copy constructor
2544OsiOldLinkBranchingObject::OsiOldLinkBranchingObject ( const OsiOldLinkBranchingObject & rhs) :OsiSOSBranchingObject(rhs)
2545{
2546}
2547
2548// Assignment operator
2549OsiOldLinkBranchingObject & 
2550OsiOldLinkBranchingObject::operator=( const OsiOldLinkBranchingObject& rhs)
2551{
2552  if (this != &rhs) {
2553    OsiSOSBranchingObject::operator=(rhs);
2554  }
2555  return *this;
2556}
2557OsiBranchingObject * 
2558OsiOldLinkBranchingObject::clone() const
2559{ 
2560  return (new OsiOldLinkBranchingObject(*this));
2561}
2562
2563
2564// Destructor
2565OsiOldLinkBranchingObject::~OsiOldLinkBranchingObject ()
2566{
2567}
2568double
2569OsiOldLinkBranchingObject::branch(OsiSolverInterface * solver)
2570{
2571  const OsiOldLink * set =
2572    dynamic_cast <const OsiOldLink *>(originalObject_) ;
2573  assert (set);
2574  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
2575  branchIndex_++;
2576  int numberMembers = set->numberMembers();
2577  const int * which = set->members();
2578  const double * weights = set->weights();
2579  int numberLinks = set->numberLinks();
2580  //const double * lower = info->lower_;
2581  //const double * upper = solver->getColUpper();
2582  // *** for way - up means fix all those in down section
2583  if (way<0) {
2584    int i;
2585    for ( i=0;i<numberMembers;i++) {
2586      if (weights[i] > value_)
2587        break;
2588    }
2589    assert (i<numberMembers);
2590    int base=i*numberLinks;;
2591    for (;i<numberMembers;i++) {
2592      for (int k=0;k<numberLinks;k++) {
2593        int iColumn = which[base+k];
2594        solver->setColUpper(iColumn,0.0);
2595      }
2596      base += numberLinks;
2597    }
2598  } else {
2599    int i;
2600    int base=0;
2601    for ( i=0;i<numberMembers;i++) { 
2602      if (weights[i] >= value_) {
2603        break;
2604      } else {
2605        for (int k=0;k<numberLinks;k++) {
2606          int iColumn = which[base+k];
2607          solver->setColUpper(iColumn,0.0);
2608        }
2609        base += numberLinks;
2610      }
2611    }
2612    assert (i<numberMembers);
2613  }
2614  return 0.0;
2615}
2616// Print what would happen 
2617void
2618OsiOldLinkBranchingObject::print(const OsiSolverInterface * solver)
2619{
2620  const OsiOldLink * set =
2621    dynamic_cast <const OsiOldLink *>(originalObject_) ;
2622  assert (set);
2623  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
2624  int numberMembers = set->numberMembers();
2625  int numberLinks = set->numberLinks();
2626  const double * weights = set->weights();
2627  const int * which = set->members();
2628  const double * upper = solver->getColUpper();
2629  int first=numberMembers;
2630  int last=-1;
2631  int numberFixed=0;
2632  int numberOther=0;
2633  int i;
2634  int base=0;
2635  for ( i=0;i<numberMembers;i++) {
2636    for (int k=0;k<numberLinks;k++) {
2637      int iColumn = which[base+k];
2638      double bound = upper[iColumn];
2639      if (bound) {
2640        first = CoinMin(first,i);
2641        last = CoinMax(last,i);
2642      }
2643    }
2644    base += numberLinks;
2645  }
2646  // *** for way - up means fix all those in down section
2647  base=0;
2648  if (way<0) {
2649    printf("SOS Down");
2650    for ( i=0;i<numberMembers;i++) {
2651      if (weights[i] > value_) 
2652        break;
2653      for (int k=0;k<numberLinks;k++) {
2654        int iColumn = which[base+k];
2655        double bound = upper[iColumn];
2656        if (bound)
2657          numberOther++;
2658      }
2659      base += numberLinks;
2660    }
2661    assert (i<numberMembers);
2662    for (;i<numberMembers;i++) {
2663      for (int k=0;k<numberLinks;k++) {
2664        int iColumn = which[base+k];
2665        double bound = upper[iColumn];
2666        if (bound)
2667          numberFixed++;
2668      }
2669      base += numberLinks;
2670    }
2671  } else {
2672    printf("SOS Up");
2673    for ( i=0;i<numberMembers;i++) {
2674      if (weights[i] >= value_)
2675        break;
2676      for (int k=0;k<numberLinks;k++) {
2677        int iColumn = which[base+k];
2678        double bound = upper[iColumn];
2679        if (bound)
2680          numberFixed++;
2681      }
2682      base += numberLinks;
2683    }
2684    assert (i<numberMembers);
2685    for (;i<numberMembers;i++) {
2686      for (int k=0;k<numberLinks;k++) {
2687        int iColumn = which[base+k];
2688        double bound = upper[iColumn];
2689        if (bound)
2690          numberOther++;
2691      }
2692      base += numberLinks;
2693    }
2694  }
2695  assert ((numberFixed%numberLinks)==0);
2696  assert ((numberOther%numberLinks)==0);
2697  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
2698         value_,first,weights[first],last,weights[last],numberFixed/numberLinks,
2699         numberOther/numberLinks);
2700}
2701// Default Constructor
2702OsiBiLinear::OsiBiLinear ()
2703  : OsiObject2(),
2704    coefficient_(0.0),
2705    xMeshSize_(0.0),
2706    yMeshSize_(0.0),
2707    xSatisfied_(1.0e-6),
2708    ySatisfied_(1.0e-6),
2709    xySatisfied_(1.0e-6),
2710    xyBranchValue_(0.0),
2711    xColumn_(-1),
2712    yColumn_(-1),
2713    firstLambda_(-1),
2714    branchingStrategy_(0),
2715    boundType_(0),
2716    xRow_(-1),
2717    yRow_(-1),
2718    xyRow_(-1),
2719    convexity_(-1),
2720    chosen_(-1)
2721{
2722}
2723
2724// Useful constructor
2725OsiBiLinear::OsiBiLinear (OsiSolverInterface * solver, int xColumn,
2726                          int yColumn, int xyRow, double coefficient,
2727                          double xMesh, double yMesh,
2728                          int numberExistingObjects,const OsiObject ** objects )
2729  : OsiObject2(),
2730    coefficient_(coefficient),
2731    xMeshSize_(xMesh),
2732    yMeshSize_(yMesh),
2733    xSatisfied_(1.0e-6),
2734    ySatisfied_(1.0e-6),
2735    xySatisfied_(1.0e-6),
2736    xyBranchValue_(0.0),
2737    xColumn_(xColumn),
2738    yColumn_(yColumn),
2739    firstLambda_(-1),
2740    branchingStrategy_(0),
2741    boundType_(0),
2742    xRow_(-1),
2743    yRow_(-1),
2744    xyRow_(xyRow),
2745    convexity_(-1),
2746    chosen_(-1)
2747{
2748  double columnLower[4];
2749  double columnUpper[4];
2750  double objective[4];
2751  double rowLower[3];
2752  double rowUpper[3];
2753  CoinBigIndex starts[5];
2754  int index[16];
2755  double element[16];
2756  int i;
2757  starts[0]=0;
2758  // rows
2759  int numberRows = solver->getNumRows();
2760  // convexity
2761  rowLower[0]=1.0;
2762  rowUpper[0]=1.0;
2763  convexity_ = numberRows;
2764  starts[1]=0;
2765  // x
2766  rowLower[1]=0.0;
2767  rowUpper[1]=0.0;
2768  index[0]=xColumn_;
2769  element[0]=-1.0;
2770  xRow_ = numberRows+1;
2771  starts[2]=1;
2772  int nAdd=2;
2773  if (xColumn_!=yColumn_) {
2774    rowLower[2]=0.0;
2775    rowUpper[2]=0.0;
2776    index[1]=yColumn;
2777    element[1]=-1.0;
2778    nAdd=3;
2779    yRow_ = numberRows+2;
2780    starts[3]=2;
2781  } else {
2782    yRow_=-1;
2783    branchingStrategy_=1;
2784  }
2785  // may be objective
2786  assert (xyRow_>=-1);
2787  solver->addRows(nAdd,starts,index,element,rowLower,rowUpper);
2788  int n=0;
2789  // order is LxLy, LxUy, UxLy and UxUy
2790  firstLambda_ = solver->getNumCols();
2791  // bit sloppy as theoretically could be infeasible but otherwise need to do more work
2792  double xB[2];
2793  double yB[2];
2794  const double * lower = solver->getColLower();
2795  const double * upper = solver->getColUpper();
2796  xB[0]=lower[xColumn_];
2797  xB[1]=upper[xColumn_];
2798  yB[0]=lower[yColumn_];
2799  yB[1]=upper[yColumn_];
2800  if (xMeshSize_!=floor(xMeshSize_)) {
2801    // not integral
2802    xSatisfied_ = CoinMax(xSatisfied_,0.51*xMeshSize_);
2803    if (!yMeshSize_) {
2804      xySatisfied_ = CoinMax(xySatisfied_,xSatisfied_*CoinMax(fabs(yB[0]),fabs(yB[1])));
2805    }
2806  }
2807  if (yMeshSize_!=floor(yMeshSize_)) {
2808    // not integral
2809    ySatisfied_ = CoinMax(ySatisfied_,0.51*yMeshSize_);
2810    if (!xMeshSize_) {
2811      xySatisfied_ = CoinMax(xySatisfied_,ySatisfied_*CoinMax(fabs(xB[0]),fabs(xB[1])));
2812    }
2813  }
2814  // adjust
2815  double distance;
2816  double steps;
2817  if (xMeshSize_) {
2818    distance = xB[1]-xB[0];
2819    steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
2820    distance = xB[0]+xMeshSize_*steps;
2821    if (fabs(xB[1]-distance)>xSatisfied_) {
2822      printf("bad x mesh %g %g %g -> %g\n",xB[0],xMeshSize_,xB[1],distance);
2823      //double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_);
2824      //printf("xSatisfied increased to %g\n",newValue);
2825      //xSatisfied_ = newValue;
2826      //xB[1]=distance;
2827      //solver->setColUpper(xColumn_,distance);
2828    }
2829  }
2830  if (yMeshSize_) {
2831    distance = yB[1]-yB[0];
2832    steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
2833    distance = yB[0]+yMeshSize_*steps;
2834    if (fabs(yB[1]-distance)>ySatisfied_) {
2835      printf("bad y mesh %g %g %g -> %g\n",yB[0],yMeshSize_,yB[1],distance);
2836      //double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_);
2837      //printf("ySatisfied increased to %g\n",newValue);
2838      //ySatisfied_ = newValue;
2839      //yB[1]=distance;
2840      //solver->setColUpper(yColumn_,distance);
2841    }
2842  }
2843  for (i=0;i<4;i++) {
2844    double x = (i<2) ? xB[0] : xB[1];
2845    double y = ((i&1)==0) ? yB[0] : yB[1];
2846    columnLower[i]=0.0;
2847    columnUpper[i]=2.0;
2848    objective[i]=0.0;
2849    double value;
2850    // xy
2851    value=coefficient_*x*y;
2852    if (xyRow_>=0) { 
2853      if (fabs(value)<1.0e-19)
2854        value = 1.0e-19;
2855      element[n]=value;
2856      index[n++]=xyRow_;
2857    } else {
2858      objective[i]=value;
2859    }
2860    // convexity
2861    value=1.0;
2862    element[n]=value;
2863    index[n++]=0+numberRows;
2864    // x
2865    value=x;
2866    if (fabs(value)<1.0e-19)
2867      value = 1.0e-19;
2868    element[n]=value;
2869    index[n++]=1+numberRows;
2870    if (xColumn_!=yColumn_) {
2871      // y
2872      value=y;
2873      if (fabs(value)<1.0e-19)
2874      value = 1.0e-19;
2875      element[n]=value;
2876      index[n++]=2+numberRows;
2877    }
2878    starts[i+1]=n;
2879  }
2880  solver->addCols(4,starts,index,element,columnLower,columnUpper,objective);
2881  // At least one has to have a mesh
2882  if (!xMeshSize_&&(!yMeshSize_||yRow_<0)) {
2883    printf("one of x and y must have a mesh size\n");
2884    abort();
2885  } else if (yRow_>=0) {
2886    if (!xMeshSize_)
2887      branchingStrategy_ = 2;
2888    else if (!yMeshSize_)
2889      branchingStrategy_ = 1;
2890  }
2891  // Now add constraints to link in x and or y to existing ones.
2892  bool xDone=false;
2893  bool yDone=false;
2894  // order is LxLy, LxUy, UxLy and UxUy
2895  for (i=numberExistingObjects-1;i>=0;i--) {
2896    const OsiObject * obj = objects[i];
2897    const OsiBiLinear * obj2 =
2898      dynamic_cast <const OsiBiLinear *>(obj) ;
2899    if (obj2) {
2900      if (xColumn_==obj2->xColumn_&&!xDone) {
2901        // make sure y equal
2902        double rhs=0.0;
2903        CoinBigIndex starts[2];
2904        int index[4];
2905        double element[4]= {1.0,1.0,-1.0,-1.0};
2906        starts[0]=0;
2907        starts[1]=4;
2908        index[0]=firstLambda_+0;
2909        index[1]=firstLambda_+1;
2910        index[2]=obj2->firstLambda_+0;
2911        index[3]=obj2->firstLambda_+1;
2912        solver->addRows(1,starts,index,element,&rhs,&rhs);
2913        xDone=true;
2914      }
2915      if (yColumn_==obj2->yColumn_&&yRow_>=0&&!yDone) {
2916        // make sure x equal
2917        double rhs=0.0;
2918        CoinBigIndex starts[2];
2919        int index[4];
2920        double element[4]= {1.0,1.0,-1.0,-1.0};
2921        starts[0]=0;
2922        starts[1]=4;
2923        index[0]=firstLambda_+0;
2924        index[1]=firstLambda_+2;
2925        index[2]=obj2->firstLambda_+0;
2926        index[3]=obj2->firstLambda_+2;
2927        solver->addRows(1,starts,index,element,&rhs,&rhs);
2928        yDone=true;
2929      }
2930    }
2931  }
2932}
2933
2934// Copy constructor
2935OsiBiLinear::OsiBiLinear ( const OsiBiLinear & rhs)
2936  :OsiObject2(rhs),
2937   coefficient_(rhs.coefficient_),
2938   xMeshSize_(rhs.xMeshSize_),
2939   yMeshSize_(rhs.yMeshSize_),
2940   xSatisfied_(rhs.xSatisfied_),
2941   ySatisfied_(rhs.ySatisfied_),
2942   xySatisfied_(rhs.xySatisfied_),
2943   xyBranchValue_(rhs.xyBranchValue_),
2944   xColumn_(rhs.xColumn_),
2945   yColumn_(rhs.yColumn_),
2946   firstLambda_(rhs.firstLambda_),
2947   branchingStrategy_(rhs.branchingStrategy_),
2948   boundType_(rhs.boundType_),
2949   xRow_(rhs.xRow_),
2950   yRow_(rhs.yRow_),
2951   xyRow_(rhs.xyRow_),
2952   convexity_(rhs.convexity_),
2953   chosen_(rhs.chosen_)
2954{
2955}
2956
2957// Clone
2958OsiObject *
2959OsiBiLinear::clone() const
2960{
2961  return new OsiBiLinear(*this);
2962}
2963
2964// Assignment operator
2965OsiBiLinear & 
2966OsiBiLinear::operator=( const OsiBiLinear& rhs)
2967{
2968  if (this!=&rhs) {
2969    OsiObject2::operator=(rhs);
2970    coefficient_ = rhs.coefficient_;
2971    xMeshSize_ = rhs.xMeshSize_;
2972    yMeshSize_ = rhs.yMeshSize_;
2973    xSatisfied_ = rhs.xSatisfied_;
2974    ySatisfied_ = rhs.ySatisfied_;
2975    xySatisfied_ = rhs.xySatisfied_;
2976    xyBranchValue_ = rhs.xyBranchValue_;
2977    xColumn_ = rhs.xColumn_;
2978    yColumn_ = rhs.yColumn_;
2979    firstLambda_ = rhs.firstLambda_;
2980    branchingStrategy_ = rhs.branchingStrategy_;
2981    boundType_ = rhs.boundType_;
2982    xRow_ = rhs.xRow_;
2983    yRow_ = rhs.yRow_;
2984    xyRow_ = rhs.xyRow_;
2985    convexity_ = rhs.convexity_;
2986    chosen_ = rhs.chosen_;
2987  }
2988  return *this;
2989}
2990
2991// Destructor
2992OsiBiLinear::~OsiBiLinear ()
2993{
2994}
2995
2996// Infeasibility - large is 0.5
2997double 
2998OsiBiLinear::infeasibility(const OsiBranchingInformation * info,int & whichWay) const
2999{
3000  // order is LxLy, LxUy, UxLy and UxUy
3001  double xB[2];
3002  double yB[2];
3003  xB[0]=info->lower_[xColumn_];
3004  xB[1]=info->upper_[xColumn_];
3005  yB[0]=info->lower_[yColumn_];
3006  yB[1]=info->upper_[yColumn_];
3007  double x = info->solution_[xColumn_];
3008  x = CoinMax(x,xB[0]);
3009  x = CoinMin(x,xB[1]);
3010  double y = info->solution_[yColumn_];
3011  y = CoinMax(y,yB[0]);
3012  y = CoinMin(y,yB[1]);
3013  int j;
3014#ifndef NDEBUG
3015  double xLambda = 0.0;
3016  double yLambda = 0.0;
3017  if ((branchingStrategy_&4)==0) {
3018    for (j=0;j<4;j++) {
3019      int iX = j>>1;
3020      int iY = j&1;
3021      xLambda += xB[iX]*info->solution_[firstLambda_+j];
3022      if (yRow_>=0)
3023        yLambda += yB[iY]*info->solution_[firstLambda_+j];
3024    }
3025  } else {
3026    const double * element = info->elementByColumn_;
3027    const int * row = info->row_;
3028    const CoinBigIndex * columnStart = info->columnStart_;
3029    const int * columnLength = info->columnLength_;
3030    for (j=0;j<4;j++) {
3031      int iColumn = firstLambda_+j;
3032      int iStart = columnStart[iColumn];
3033      int iEnd = iStart + columnLength[iColumn];
3034      int k=iStart;
3035      double sol = info->solution_[iColumn];
3036      for (;k<iEnd;k++) {
3037        if (xRow_==row[k])
3038          xLambda += element[k]*sol;
3039        if (yRow_==row[k])
3040          yLambda += element[k]*sol;
3041      }
3042    }
3043  }
3044  assert (fabs(x-xLambda)<1.0e-1);
3045  if (yRow_>=0)
3046    assert (fabs(y-yLambda)<1.0e-1);
3047#endif
3048  // If x or y not satisfied then branch on that
3049  double distance;
3050  double steps;
3051  bool xSatisfied;
3052  double xNew;
3053  if (xMeshSize_) {
3054    if (x<0.5*(xB[0]+xB[1])) {
3055      distance = x-xB[0];
3056      steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
3057      xNew = xB[0]+steps*xMeshSize_;
3058      assert (xNew<=xB[1]+xSatisfied_);
3059      xSatisfied =  (fabs(xNew-x)<xSatisfied_);
3060    } else {
3061      distance = xB[1]-x;
3062      steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
3063      xNew = xB[1]-steps*xMeshSize_;
3064      assert (xNew>=xB[0]-xSatisfied_);
3065      xSatisfied =  (fabs(xNew-x)<xSatisfied_);
3066    }
3067  } else {
3068    xSatisfied=true;
3069  }
3070  bool ySatisfied;
3071  double yNew;
3072  if (yMeshSize_) {
3073    if (y<0.5*(yB[0]+yB[1])) {
3074      distance = y-yB[0];
3075      steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
3076      yNew = yB[0]+steps*yMeshSize_;
3077      assert (yNew<=yB[1]+ySatisfied_);
3078      ySatisfied =  (fabs(yNew-y)<ySatisfied_);
3079    } else {
3080      distance = yB[1]-y;
3081      steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
3082      yNew = yB[1]-steps*yMeshSize_;
3083      assert (yNew>=yB[0]-ySatisfied_);
3084      ySatisfied =  (fabs(yNew-y)<ySatisfied_);
3085    }
3086  } else {
3087    ySatisfied=true;
3088  }
3089  /* There are several possibilities
3090     1 - one or both are unsatisfied and branching strategy tells us what to do
3091     2 - both are unsatisfied and branching strategy is 0
3092     3 - both are satisfied but xy is not
3093         3a one has bounds within satisfied_ - other does not
3094         (or neither have but branching strategy tells us what to do)
3095         3b neither do - and branching strategy does not tell us
3096         3c both do - treat as feasible knowing another copy of object will fix
3097     4 - both are satisfied and xy is satisfied - as 3c
3098  */
3099  chosen_=-1;
3100  xyBranchValue_=COIN_DBL_MAX;
3101  whichWay_=0;
3102  double xyTrue = x*y;
3103  double xyLambda = 0.0;
3104  if ((branchingStrategy_&4)==0) {
3105    for (j=0;j<4;j++) {
3106      int iX = j>>1;
3107      int iY = j&1;
3108      xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j];
3109    }
3110  } else {
3111    if (xyRow_>=0) {
3112      const double * element = info->elementByColumn_;
3113      const int * row = info->row_;
3114      const CoinBigIndex * columnStart = info->columnStart_;
3115      const int * columnLength = info->columnLength_;
3116      for (j=0;j<4;j++) {
3117        int iColumn = firstLambda_+j;
3118        int iStart = columnStart[iColumn];
3119        int iEnd = iStart + columnLength[iColumn];
3120        int k=iStart;
3121        double sol = info->solution_[iColumn];
3122        for (;k<iEnd;k++) {
3123          if (xyRow_==row[k])
3124            xyLambda += element[k]*sol;
3125        }
3126      }
3127    } else {
3128      // objective
3129      const double * objective = info->objective_;
3130      for (j=0;j<4;j++) {
3131        int iColumn = firstLambda_+j;
3132        double sol = info->solution_[iColumn];
3133        xyLambda += objective[iColumn]*sol;
3134      }
3135    }
3136    xyLambda /= coefficient_;
3137  }
3138  // If pseudo shadow prices then see what would happen
3139  //double pseudoEstimate = 0.0;
3140  if (info->defaultDual_>=0.0) {
3141    // If we move to xy then we move by coefficient * (xyTrue-xyLambda) on row xyRow_
3142    double move = xyTrue-xyLambda;
3143    assert (xyRow_>=0);
3144    if (boundType_==0) {
3145      move *= coefficient_;
3146      move *= info->pi_[xyRow_];
3147      move = CoinMax(move,0.0);
3148    } else if (boundType_==1) {
3149      // if OK then say satisfied
3150    } else if (boundType_==2) {
3151    } else {
3152      // == row so move x and y not xy
3153    }
3154  }
3155  if ( !xSatisfied) {
3156    if (!ySatisfied) {
3157      if ((branchingStrategy_&3)==0) {
3158        // If pseudo shadow prices then see what would happen
3159        if (info->defaultDual_>=0.0) {
3160          // need coding here
3161          if (fabs(x-xNew)>fabs(y-yNew)) {
3162            chosen_=0;
3163            xyBranchValue_=x;
3164          } else {
3165            chosen_=1;
3166            xyBranchValue_=y;
3167          }
3168        } else {
3169          if (fabs(x-xNew)>fabs(y-yNew)) {
3170            chosen_=0;
3171            xyBranchValue_=x;
3172          } else {
3173            chosen_=1;
3174            xyBranchValue_=y;
3175          }
3176        }
3177      } else if ((branchingStrategy_&3)==1) {
3178        chosen_=0;
3179        xyBranchValue_=x;
3180      } else {
3181        chosen_=1;
3182        xyBranchValue_=y;
3183      }
3184    } else {
3185      // y satisfied
3186      chosen_=0;
3187      xyBranchValue_=x;
3188    }
3189  } else {
3190    // x satisfied
3191    if (!ySatisfied) {
3192      chosen_=1;
3193      xyBranchValue_=y;
3194    } else {
3195      /*
3196        3 - both are satisfied but xy is not
3197         3a one has bounds within satisfied_ - other does not
3198         (or neither have but branching strategy tells us what to do)
3199         3b neither do - and branching strategy does not tell us
3200         3c both do - treat as feasible knowing another copy of object will fix
3201        4 - both are satisfied and xy is satisfied - as 3c
3202      */
3203      if (fabs(xyLambda-xyTrue)<xySatisfied_||(xB[0]==xB[1]&&yB[0]==yB[1])) {
3204        // satisfied
3205#if 0
3206        printf("all satisfied true %g lambda %g\n",
3207               xyTrue,xyLambda);
3208        printf("x %d (%g,%g,%g) y %d (%g,%g,%g)\n",
3209               xColumn_,xB[0],x,xB[1],
3210               yColumn_,yB[0],y,yB[1]);
3211#endif
3212      } else {
3213        // May be infeasible - check
3214        bool feasible=true;
3215        if (xB[0]==xB[1]&&yB[0]==yB[1]) {
3216          double lambda[4];
3217          computeLambdas(info->solver_, lambda);
3218          for (int j=0;j<4;j++) {
3219            int iColumn = firstLambda_+j;
3220            if (info->lower_[iColumn]>lambda[j]+1.0e-5||
3221                info->upper_[iColumn]<lambda[j]-1.0e-5)
3222              feasible=false;
3223          }
3224        }
3225        if (feasible) {
3226          if (xB[1]-xB[0]>=xSatisfied_&&xMeshSize_) {
3227            if (yB[1]-yB[0]>=ySatisfied_&&yMeshSize_) {
3228              if ((branchingStrategy_&3)==0) {
3229                // If pseudo shadow prices then see what would happen
3230                if (info->defaultDual_>=0.0) {
3231                  // need coding here
3232                  if (xB[1]-xB[0]>yB[1]-yB[0]) {
3233                    chosen_=0;
3234                    xyBranchValue_=0.5*(xB[0]+xB[1]);
3235                  } else {
3236                    chosen_=1;
3237                    xyBranchValue_=0.5*(yB[0]+yB[1]);
3238                  }
3239                } else {
3240                  if (xB[1]-xB[0]>yB[1]-yB[0]) {
3241                    chosen_=0;
3242                    xyBranchValue_=0.5*(xB[0]+xB[1]);
3243                  } else {
3244                    chosen_=1;
3245                    xyBranchValue_=0.5*(yB[0]+yB[1]);
3246                  }
3247                }
3248              } else if ((branchingStrategy_&3)==1) {
3249                chosen_=0;
3250                xyBranchValue_=0.5*(xB[0]+xB[1]);
3251              } else {
3252                chosen_=1;
3253                xyBranchValue_=0.5*(yB[0]+yB[1]);
3254              }
3255            } else {
3256              // y satisfied
3257              chosen_=0;
3258              xyBranchValue_=0.5*(xB[0]+xB[1]);
3259            }
3260          } else if (yB[1]-yB[0]>=ySatisfied_&&yMeshSize_) {
3261            chosen_=1;
3262            xyBranchValue_=0.5*(yB[0]+yB[1]);
3263          } else {
3264            // treat as satisfied unless no coefficient tightening
3265            if ((branchingStrategy_&4)!=0) {
3266              chosen_=0; // fix up in branch
3267              xyBranchValue_=x;
3268            }
3269          }
3270        } else {
3271          // node not feasible!!!
3272          chosen_=0;
3273          infeasibility_=COIN_DBL_MAX;
3274          otherInfeasibility_ = COIN_DBL_MAX;
3275          whichWay=whichWay_;
3276          return infeasibility_;
3277        }
3278      }
3279    }
3280  }
3281  if (chosen_==-1) {
3282    infeasibility_=0.0;
3283  } else if (chosen_==0) {
3284    infeasibility_ = CoinMax(fabs(xyBranchValue_-x),1.0e-12);
3285    //assert (xyBranchValue_>=info->lower_[xColumn_]&&xyBranchValue_<=info->upper_[xColumn_]);
3286  } else {
3287    infeasibility_ = CoinMax(fabs(xyBranchValue_-y),1.0e-12);
3288    //assert (xyBranchValue_>=info->lower_[yColumn_]&&xyBranchValue_<=info->upper_[yColumn_]);
3289  }
3290  if (info->defaultDual_<0.0) {
3291    // not using pseudo shadow prices
3292    otherInfeasibility_ = 1.0-infeasibility_;
3293  } else {
3294    abort();
3295  }
3296  if (infeasibility_) {
3297    bool fixed=true;
3298    for (int j=0;j<4;j++) {
3299      int iColumn = firstLambda_+j;
3300      //if (info->lower_[iColumn]) printf("lower of %g on %d\n",info->lower_[iColumn],iColumn);
3301      if (info->lower_[iColumn]<info->upper_[iColumn])
3302        fixed=false;
3303    }
3304    if (fixed) {
3305      //printf("must be tolerance problem - xy true %g lambda %g\n",xyTrue,xyLambda);
3306      chosen_=-1;
3307      infeasibility_=0.0;
3308    }
3309  }
3310  whichWay=whichWay_;
3311  return infeasibility_;
3312}
3313
3314// This looks at solution and sets bounds to contain solution
3315double
3316OsiBiLinear::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
3317{
3318  // If another object has finer mesh ignore this
3319  if ((branchingStrategy_&8)!=0)
3320    return 0.0;
3321  // order is LxLy, LxUy, UxLy and UxUy
3322  double xB[2];
3323  double yB[2];
3324  xB[0]=info->lower_[xColumn_];
3325  xB[1]=info->upper_[xColumn_];
3326  yB[0]=info->lower_[yColumn_];
3327  yB[1]=info->upper_[yColumn_];
3328  double x = info->solution_[xColumn_];
3329  double y = info->solution_[yColumn_];
3330  int j;
3331#ifndef NDEBUG
3332  double xLambda = 0.0;
3333  double yLambda = 0.0;
3334  if ((branchingStrategy_&4)==0) {
3335    for (j=0;j<4;j++) {
3336      int iX = j>>1;
3337      int iY = j&1;
3338      xLambda += xB[iX]*info->solution_[firstLambda_+j];
3339      if (yRow_>=0)
3340        yLambda += yB[iY]*info->solution_[firstLambda_+j];
3341    }
3342  } else {
3343    const double * element = info->elementByColumn_;
3344    const int * row = info->row_;
3345    const CoinBigIndex * columnStart = info->columnStart_;
3346    const int * columnLength = info->columnLength_;
3347    for (j=0;j<4;j++) {
3348      int iColumn = firstLambda_+j;
3349      int iStart = columnStart[iColumn];
3350      int iEnd = iStart + columnLength[iColumn];
3351      int k=iStart;
3352      double sol = info->solution_[iColumn];
3353      for (;k<iEnd;k++) {
3354        if (xRow_==row[k])
3355          xLambda += element[k]*sol;
3356        if (yRow_==row[k])
3357          yLambda += element[k]*sol;
3358      }
3359    }
3360  }
3361  if (yRow_<0)
3362    yLambda = xLambda;
3363#if 0
3364  if (fabs(x-xLambda)>1.0e-4||
3365      fabs(y-yLambda)>1.0e-4)
3366    printf("feasibleregion x %d given %g lambda %g y %d given %g lambda %g\n",
3367           xColumn_,x,xLambda,
3368           yColumn_,y,yLambda);
3369#endif
3370#endif
3371  double infeasibility=0.0;
3372  double distance;
3373  double steps;
3374  double xNew=x;
3375  if (xMeshSize_) {
3376    distance = x-xB[0];
3377    if (x<0.5*(xB[0]+xB[1])) {
3378      distance = x-xB[0];
3379      steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
3380      xNew = xB[0]+steps*xMeshSize_;
3381      assert (xNew<=xB[1]+xSatisfied_);
3382    } else {
3383      distance = xB[1]-x;
3384      steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
3385      xNew = xB[1]-steps*xMeshSize_;
3386      assert (xNew>=xB[0]-xSatisfied_);
3387    }
3388    if (xMeshSize_<1.0&&fabs(xNew-x)<=xSatisfied_) {
3389      double lo = CoinMax(xB[0],x-0.5*xSatisfied_);
3390      double up = CoinMin(xB[1],x+0.5*xSatisfied_);
3391      solver->setColLower(xColumn_,lo);
3392      solver->setColUpper(xColumn_,up);
3393    } else {
3394      infeasibility +=  fabs(xNew-x);
3395      solver->setColLower(xColumn_,xNew);
3396      solver->setColUpper(xColumn_,xNew);
3397    }
3398  }
3399  double yNew=y;
3400  if (yMeshSize_) {
3401    distance = y-yB[0];
3402    if (y<0.5*(yB[0]+yB[1])) {
3403      distance = y-yB[0];
3404      steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
3405      yNew = yB[0]+steps*yMeshSize_;
3406      assert (yNew<=yB[1]+ySatisfied_);
3407    } else {
3408      distance = yB[1]-y;
3409      steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
3410      yNew = yB[1]-steps*yMeshSize_;
3411      assert (yNew>=yB[0]-ySatisfied_);
3412    }
3413    if (yMeshSize_<1.0&&fabs(yNew-y)<=ySatisfied_) {
3414      double lo = CoinMax(yB[0],y-0.5*ySatisfied_);
3415      double up = CoinMin(yB[1],y+0.5*ySatisfied_);
3416      solver->setColLower(yColumn_,lo);
3417      solver->setColUpper(yColumn_,up);
3418    } else {
3419      infeasibility +=  fabs(yNew-y);
3420      solver->setColLower(yColumn_,yNew);
3421      solver->setColUpper(yColumn_,yNew);
3422    }
3423  } 
3424  if (0) {
3425    // temp
3426      solver->setColLower(xColumn_,x);
3427      solver->setColUpper(xColumn_,x);
3428      solver->setColLower(yColumn_,y);
3429      solver->setColUpper(yColumn_,y);
3430  }
3431  if ((branchingStrategy_&4)) {
3432    // fake to make correct
3433    double lambda[4];
3434    computeLambdas(solver,lambda);
3435    for (int j=0;j<4;j++) {
3436      int iColumn = firstLambda_+j;
3437      double value = lambda[j];
3438      solver->setColLower(iColumn,value);
3439      solver->setColUpper(iColumn,value);
3440    }
3441  }
3442  double xyTrue = xNew*yNew;
3443  double xyLambda = 0.0;
3444  for (j=0;j<4;j++) {
3445    int iX = j>>1;
3446    int iY = j&1;
3447    xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j];
3448  }
3449  infeasibility += fabs(xyTrue-xyLambda);
3450  return infeasibility;
3451}
3452// Returns true value of single xyRow coefficient
3453double 
3454OsiBiLinear::xyCoefficient(const double * solution) const
3455{
3456  // If another object has finer mesh ignore this
3457  if ((branchingStrategy_&8)!=0)
3458    return 0.0;
3459  double x = solution[xColumn_];
3460  double y = solution[yColumn_];
3461  //printf("x (%d,%g) y (%d,%g) x*y*coefficient %g\n",
3462  // xColumn_,x,yColumn_,y,x*y*coefficient_);
3463  return x*y*coefficient_;
3464}
3465
3466// Redoes data when sequence numbers change
3467void 
3468OsiBiLinear::resetSequenceEtc(int numberColumns, const int * originalColumns)
3469{
3470  int i;
3471  for (i=0;i<numberColumns;i++) {
3472    if (originalColumns[i]==firstLambda_)
3473      break;
3474  }
3475  if (i<numberColumns) {
3476    firstLambda_ = i;
3477    for (int j=0;j<4;j++) {
3478      assert (originalColumns[j+i]-firstLambda_==j);
3479    }
3480  } else {
3481    printf("lost set\n");
3482    abort();
3483  }
3484  // rows will be out anyway
3485  abort();
3486}
3487
3488// Creates a branching object
3489OsiBranchingObject * 
3490OsiBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
3491{
3492  // create object
3493  OsiBranchingObject * branch;
3494  assert (chosen_==0||chosen_==1);
3495  //if (chosen_==0)
3496  //assert (xyBranchValue_>=info->lower_[xColumn_]&&xyBranchValue_<=info->upper_[xColumn_]);
3497  //else
3498  //assert (xyBranchValue_>=info->lower_[yColumn_]&&xyBranchValue_<=info->upper_[yColumn_]);
3499  branch = new OsiBiLinearBranchingObject(solver,this,way,xyBranchValue_,chosen_);
3500  return branch;
3501}
3502// Does work of branching
3503void 
3504OsiBiLinear::newBounds(OsiSolverInterface * solver, int way, short xOrY, double separator) const
3505{
3506  int iColumn;
3507  double mesh;
3508  double satisfied;
3509  if (xOrY==0) {
3510    iColumn=xColumn_;
3511    mesh=xMeshSize_;
3512    satisfied=xSatisfied_;
3513  } else {
3514    iColumn=yColumn_;
3515    mesh=yMeshSize_;
3516    satisfied=ySatisfied_;
3517  }
3518  const double * columnLower = solver->getColLower();
3519  const double * columnUpper = solver->getColUpper();
3520  double lower = columnLower[iColumn];
3521  double distance;
3522  double steps;
3523  double zNew=separator;
3524  distance = separator-lower;
3525  assert (mesh);
3526  steps = floor ((distance+0.5*mesh)/mesh);
3527  if (mesh<1.0)
3528    zNew = lower+steps*mesh;
3529  if (zNew>columnUpper[iColumn]-satisfied)
3530    zNew = 0.5*(columnUpper[iColumn]-lower);
3531  double oldUpper = columnUpper[iColumn] ;
3532  double oldLower = columnLower[iColumn] ;
3533#ifndef NDEBUG
3534  int nullChange=0;
3535#endif
3536  if (way<0) {
3537    if (zNew>separator&&mesh<1.0)
3538      zNew -= mesh;
3539    double oldUpper = columnUpper[iColumn] ;
3540    if (zNew+satisfied>=oldUpper)
3541      zNew =0.5*(oldUpper+oldLower);
3542    if (mesh==1.0)
3543      zNew = floor(separator);
3544#ifndef NDEBUG
3545    if (oldUpper<zNew+1.0e-8)
3546      nullChange=-1;
3547#endif
3548    solver->setColUpper(iColumn,zNew);
3549  } else {
3550    if (zNew<separator&&mesh<1.0)
3551      zNew += mesh;
3552    if (zNew-satisfied<=oldLower)
3553      zNew =0.5*(oldUpper+oldLower);
3554    if (mesh==1.0)
3555      zNew = ceil(separator);
3556#ifndef NDEBUG
3557    if (oldLower>zNew-1.0e-8)
3558      nullChange=1;
3559#endif
3560    solver->setColLower(iColumn,zNew);
3561  }
3562  if ((branchingStrategy_&4)!=0
3563      &&columnLower[xColumn_]==columnUpper[xColumn_]
3564      &&columnLower[yColumn_]==columnUpper[yColumn_]) {
3565    // fake to make correct
3566    double lambda[4];
3567    computeLambdas(solver,lambda);
3568    for (int j=0;j<4;j++) {
3569      int iColumn = firstLambda_+j;
3570      double value = lambda[j];
3571#ifndef NDEBUG
3572      if (fabs(value-columnLower[iColumn])>1.0e-5||
3573          fabs(value-columnUpper[iColumn])>1.0e-5)
3574        nullChange=0;
3575#endif
3576      solver->setColLower(iColumn,value);
3577      solver->setColUpper(iColumn,value);
3578    }
3579  }
3580#ifndef NDEBUG
3581  if (nullChange)
3582    printf("null change on column%s %d - bounds %g,%g\n",nullChange>0 ? "Lower" : "Upper",
3583           iColumn,oldLower,oldUpper);
3584#endif
3585#if 0
3586  // always free up lambda
3587  for (int i=firstLambda_;i<firstLambda_+4;i++) {
3588    solver->setColLower(i,0.0);
3589    solver->setColUpper(i,2.0);
3590  }
3591#endif
3592  double xB[3];
3593  xB[0]=columnLower[xColumn_];
3594  xB[1]=columnUpper[xColumn_];
3595  double yB[3];
3596  yB[0]=columnLower[yColumn_];
3597  yB[1]=columnUpper[yColumn_];
3598  if (false&&(branchingStrategy_&4)!=0&&yRow_>=0&&
3599      xMeshSize_==1.0&&yMeshSize_==1.0) {
3600    if ((xB[1]-xB[0])*(yB[1]-yB[0])<40) {
3601      // try looking at all solutions
3602      double lower[4];
3603      double upper[4];
3604      double lambda[4];
3605      int i;
3606      double lowerLambda[4];
3607      double upperLambda[4];
3608      for (i=0;i<4;i++) {
3609        lower[i] = CoinMax(0.0,columnLower[firstLambda_+i]);
3610        upper[i] = CoinMin(1.0,columnUpper[firstLambda_+i]);
3611        lowerLambda[i]=1.0;
3612        upperLambda[i]=0.0;
3613      }
3614      // get coefficients
3615      double xybar[4];
3616      getCoefficients(solver,xB,yB,xybar);
3617      double x,y;
3618      for (x=xB[0];x<=xB[1];x++) {
3619        xB[2]=x;
3620        for (y=yB[0];y<=yB[1];y++) {
3621          yB[2]=y;
3622          computeLambdas(xB, yB, xybar,lambda);
3623          for (i=0;i<4;i++) {
3624            lowerLambda[i] = CoinMin(lowerLambda[i],lambda[i]);
3625            upperLambda[i] = CoinMax(upperLambda[i],lambda[i]);
3626          }
3627        }
3628      }
3629      double change=0.0;;
3630      for (i=0;i<4;i++) {
3631        if (lowerLambda[i]>lower[i]+1.0e-12) {
3632          solver->setColLower(firstLambda_+i,lowerLambda[i]);
3633          change += lowerLambda[i]-lower[i];
3634        }
3635        if (upperLambda[i]<upper[i]-1.0e-12) {
3636          solver->setColUpper(firstLambda_+i,upperLambda[i]);
3637          change -= upperLambda[i]-upper[i];
3638        }
3639      }
3640      if (change>1.0e-5)
3641        printf("change of %g\n",change);
3642    }
3643  }
3644  if (boundType_) {
3645    assert (!xMeshSize_||!yMeshSize_);
3646    if (xMeshSize_) {
3647      // can tighten bounds on y
3648      if ((boundType_&1)!=0) {
3649        if (xB[0]*yB[1]>coefficient_) {
3650          // tighten upper bound on y
3651          solver->setColUpper(yColumn_,coefficient_/xB[0]);
3652        }
3653      }
3654      if ((boundType_&2)!=0) {
3655        if (xB[1]*yB[0]<coefficient_) {
3656          // tighten lower bound on y
3657          solver->setColLower(yColumn_,coefficient_/xB[1]);
3658        }
3659      }
3660    } else {
3661      // can tighten bounds on x
3662      if ((boundType_&1)!=0) {
3663        if (yB[0]*xB[1]>coefficient_) {
3664          // tighten upper bound on x
3665          solver->setColUpper(xColumn_,coefficient_/yB[0]);
3666        }
3667      }
3668      if ((boundType_&2)!=0) {
3669        if (yB[1]*xB[0]<coefficient_) {
3670          // tighten lower bound on x
3671          solver->setColLower(xColumn_,coefficient_/yB[1]);
3672        }
3673      }
3674    }
3675  }
3676}
3677// Compute lambdas if coefficients not changing
3678void 
3679OsiBiLinear::computeLambdas(const OsiSolverInterface * solver,double lambda[4]) const
3680{
3681  // fix so correct
3682  double xB[3],yB[3];
3683  double xybar[4];
3684  getCoefficients(solver,xB,yB,xybar);
3685  double x,y;
3686  x=solver->getColLower()[xColumn_];
3687  assert(x==solver->getColUpper()[xColumn_]);
3688  xB[2]=x;
3689  y=solver->getColLower()[yColumn_];
3690  assert(y==solver->getColUpper()[yColumn_]);
3691  yB[2]=y;
3692  computeLambdas(xB, yB, xybar,lambda);
3693  assert (xyRow_>=0);
3694}
3695// Get LU coefficients from matrix
3696void 
3697OsiBiLinear::getCoefficients(const OsiSolverInterface * solver,double xB[2], double yB[2],
3698                             double xybar[4]) const
3699{
3700  const CoinPackedMatrix * matrix = solver->getMatrixByCol();
3701  const double * element = matrix->getElements();
3702  const double * objective = solver->getObjCoefficients();
3703  const int * row = matrix->getIndices();
3704  const CoinBigIndex * columnStart = matrix->getVectorStarts();
3705  const int * columnLength = matrix->getVectorLengths();
3706  // order is LxLy, LxUy, UxLy and UxUy
3707  int j;
3708  double multiplier = (boundType_==0) ? 1.0/coefficient_ : 1.0;
3709  if (yRow_>=0) {
3710    for (j=0;j<4;j++) {
3711      int iColumn = firstLambda_+j;
3712      int iStart = columnStart[iColumn];
3713      int iEnd = iStart + columnLength[iColumn];
3714      int k=iStart;
3715      double x=0.0;
3716      double y=0.0;
3717      xybar[j]=0.0;
3718      for (;k<iEnd;k++) {
3719        if (xRow_==row[k])
3720          x = element[k];
3721        if (yRow_==row[k])
3722          y = element[k];
3723        if (xyRow_==row[k])
3724          xybar[j] = element[k]*multiplier;
3725      }
3726      if (xyRow_<0)
3727        xybar[j] = objective[iColumn]*multiplier;
3728      if (j==0)
3729        xB[0]=x;
3730      else if (j==1)
3731        yB[1]=y;
3732      else if (j==2)
3733        yB[0]=y;
3734      else if (j==3)
3735        xB[1]=x;
3736      assert (fabs(xybar[j]-x*y)<1.0e-4);
3737    }
3738  } else {
3739    // x==y
3740    for (j=0;j<4;j++) {
3741      int iColumn = firstLambda_+j;
3742      int iStart = columnStart[iColumn];
3743      int iEnd = iStart + columnLength[iColumn];
3744      int k=iStart;
3745      double x=0.0;
3746      xybar[j]=0.0;
3747      for (;k<iEnd;k++) {
3748        if (xRow_==row[k])
3749          x = element[k];
3750        if (xyRow_==row[k])
3751          xybar[j] = element[k]*multiplier;
3752      }
3753      if (xyRow_<0)
3754        xybar[j] = objective[iColumn]*multiplier;
3755      if (j==0) {
3756        xB[0]=x;
3757        yB[0]=x;
3758      } else if (j==2) {
3759        xB[1]=x;
3760        yB[1]=x;
3761      }
3762    }
3763    assert (fabs(xybar[0]-xB[0]*yB[0])<1.0e-4);
3764    assert (fabs(xybar[1]-xB[0]*yB[1])<1.0e-4);
3765    assert (fabs(xybar[2]-xB[1]*yB[0])<1.0e-4);
3766    assert (fabs(xybar[3]-xB[1]*yB[1])<1.0e-4);
3767  }
3768}
3769// Compute lambdas (third entry in each .B is current value)
3770double
3771OsiBiLinear::computeLambdas(const double xB[3], const double yB[3], const double xybar[4], double lambda[4]) const
3772{
3773  // fake to make correct
3774  double x = xB[2];
3775  double y = yB[2];
3776  // order is LxLy, LxUy, UxLy and UxUy
3777  // l0 + l1 = this
3778  double rhs1 = (xB[1]-x)/(xB[1]-xB[0]);
3779  // l0 + l2 = this
3780  double rhs2 = (yB[1]-y)/(yB[1]-yB[0]);
3781  // For xy (taking out l3)
3782  double rhs3 = xB[1]*yB[1]-x*y;
3783  double a0 = xB[1]*yB[1] - xB[0]*yB[0];
3784  double a1 = xB[1]*yB[1] - xB[0]*yB[1];
3785  double a2 = xB[1]*yB[1] - xB[1]*yB[0];
3786  // divide through to get l0 coefficient
3787  rhs3 /= a0;
3788  a1 /= a0;
3789  a2 /= a0;
3790  // subtract out l0
3791  double b[2][2];
3792  double rhs[2];
3793  // first for l1 and l2
3794  b[0][0] = 1.0-a1;
3795  b[0][1] = -a2;
3796  rhs[0] = rhs1 - rhs3;
3797  // second
3798  b[1][0] = -a1;
3799  b[1][1] = 1.0-a2;
3800  rhs[1] = rhs2 - rhs3;
3801  if (fabs(b[0][0])>fabs(b[0][1])) {
3802    double sub = b[1][0]/b[0][0];
3803    b[1][1] -= sub*b[0][1];
3804    rhs[1] -= sub*rhs[0];
3805    assert (fabs(b[1][1])>1.0e-12);
3806    lambda[2] = rhs[1]/b[1][1];
3807    lambda[0] = rhs2 - lambda[2];
3808    lambda[1] = rhs1 - lambda[0];
3809  } else {
3810    double sub = b[1][1]/b[0][1];
3811    b[1][0] -= sub*b[0][0];
3812    rhs[1] -= sub*rhs[0];
3813    assert (fabs(b[1][0])>1.0e-12);
3814    lambda[1] = rhs[1]/b[1][0];
3815    lambda[0] = rhs1 - lambda[1];
3816    lambda[2] = rhs2 - lambda[0];
3817  }
3818  lambda[3] = 1.0- (lambda[0]+lambda[1]+lambda[2]);
3819  double infeasibility=0.0;
3820  double xy=0.0;
3821  for (int j=0;j<4;j++) {
3822    double value = lambda[j];
3823    if (value>1.0) {
3824      infeasibility += value-1.0;
3825      value=1.0;
3826    }
3827    if (value<0.0) {
3828      infeasibility -= value;
3829      value=0.0;
3830    }
3831    lambda[j]=value;
3832    xy += xybar[j]*value;
3833  }
3834  assert (fabs(xy-x*y)<1.0e-4);
3835  return infeasibility;
3836}
3837// Updates coefficients
3838int
3839OsiBiLinear::updateCoefficients(const double * lower, const double * upper, double * objective, 
3840                                CoinPackedMatrix * matrix, CoinWarmStartBasis * basis) const
3841{
3842  // Return if no updates
3843  if ((branchingStrategy_&4)!=0)
3844    return 0;
3845  int numberUpdated=0;
3846  double * element = matrix->getMutableElements();
3847  const int * row = matrix->getIndices();
3848  const CoinBigIndex * columnStart = matrix->getVectorStarts();
3849  //const int * columnLength = matrix->getVectorLengths();
3850  // order is LxLy, LxUy, UxLy and UxUy
3851  double xB[2];
3852  double yB[2];
3853  xB[0]=lower[xColumn_];
3854  xB[1]=upper[xColumn_];
3855  yB[0]=lower[yColumn_];
3856  yB[1]=upper[yColumn_];
3857  //printf("x %d (%g,%g) y %d (%g,%g)\n",
3858  // xColumn_,xB[0],xB[1],
3859  // yColumn_,yB[0],yB[1]);
3860  CoinWarmStartBasis::Status status[4];
3861  int numStruct = basis ? basis->getNumStructural()-firstLambda_ : 0;
3862  double coefficient = (boundType_==0) ? coefficient_ : 1.0;
3863  for (int j=0;j<4;j++) {
3864    status[j]=(j<numStruct) ? basis->getStructStatus(j+firstLambda_) : CoinWarmStartBasis::atLowerBound;
3865    int iX = j>>1;
3866    double x = xB[iX];
3867    int iY = j&1;
3868    double y = yB[iY];
3869    CoinBigIndex k = columnStart[j+firstLambda_];
3870    double value;
3871    // xy
3872    value=coefficient*x*y;
3873    if (xyRow_>=0) {
3874      assert (row[k]==xyRow_);
3875#if BI_PRINT > 1
3876      printf("j %d xy (%d,%d) coeff from %g to %g\n",j,xColumn_,yColumn_,element[k],value);
3877#endif
3878      element[k++]=value;
3879    } else {
3880      // objective
3881      objective[j+firstLambda_]=value;
3882    }
3883    numberUpdated++;
3884    // convexity
3885    assert (row[k]==convexity_);
3886    k++;
3887    // x
3888    value=x;
3889#if BI_PRINT > 1
3890    printf("j %d x (%d) coeff from %g to %g\n",j,xColumn_,element[k],value);
3891#endif
3892    assert (row[k]==xRow_);
3893    element[k++]=value;
3894    numberUpdated++;
3895    if (yRow_>=0) {
3896      // y
3897      value=y;
3898#if BI_PRINT > 1
3899      printf("j %d y (%d) coeff from %g to %g\n",j,yColumn_,element[k],value);
3900#endif
3901      assert (row[k]==yRow_);
3902      element[k++]=value;
3903      numberUpdated++;
3904    }
3905  }
3906 
3907  if (xB[0]==xB[1]) {
3908    if (yB[0]==yB[1]) {
3909      // only one basic
3910      bool first=true;
3911      for (int j=0;j<4;j++) {
3912        if (status[j]==CoinWarmStartBasis::basic) {
3913          if (first) {
3914            first=false;
3915          } else {
3916            basis->setStructStatus(j+firstLambda_,CoinWarmStartBasis::atLowerBound);
3917#if BI_PRINT
3918            printf("zapping %d (x=%d,y=%d)\n",j,xColumn_,yColumn_);
3919#endif
3920          }
3921        }
3922      }
3923    } else {
3924      if (status[0]==CoinWarmStartBasis::basic&&
3925          status[2]==CoinWarmStartBasis::basic) {
3926        basis->setStructStatus(2+firstLambda_,CoinWarmStartBasis::atLowerBound);
3927#if BI_PRINT
3928        printf("zapping %d (x=%d,y=%d)\n",2,xColumn_,yColumn_);
3929#endif
3930      }
3931      if (status[1]==CoinWarmStartBasis::basic&&
3932          status[3]==CoinWarmStartBasis::basic) {
3933        basis->setStructStatus(3+firstLambda_,CoinWarmStartBasis::atLowerBound);
3934#if BI_PRINT
3935        printf("zapping %d (x=%d,y=%d)\n",3,xColumn_,yColumn_);
3936#endif
3937      }
3938    }
3939  } else if (yB[0]==yB[1]) {
3940    if (status[0]==CoinWarmStartBasis::basic&&
3941        status[1]==CoinWarmStartBasis::basic) {
3942      basis->setStructStatus(1+firstLambda_,CoinWarmStartBasis::atLowerBound);
3943#if BI_PRINT
3944      printf("zapping %d (x=%d,y=%d)\n",1,xColumn_,yColumn_);
3945#endif
3946    }
3947    if (status[2]==CoinWarmStartBasis::basic&&
3948        status[3]==CoinWarmStartBasis::basic) {
3949      basis->setStructStatus(3+firstLambda_,CoinWarmStartBasis::atLowerBound);
3950#if BI_PRINT
3951      printf("zapping %d (x=%d,y=%d)\n",3,xColumn_,yColumn_);
3952#endif
3953    }
3954  }
3955  return numberUpdated;
3956}
3957// This does NOT set mutable stuff
3958double 
3959OsiBiLinear::checkInfeasibility(const OsiBranchingInformation * info) const
3960{
3961  // If another object has finer mesh ignore this
3962  if ((branchingStrategy_&8)!=0)
3963    return 0.0;
3964  int way;
3965  double saveInfeasibility = infeasibility_;
3966  int saveWhichWay = whichWay_;
3967  double saveXyBranchValue = xyBranchValue_;
3968  short saveChosen = chosen_;
3969  double value = infeasibility(info,way);
3970  infeasibility_ = saveInfeasibility;
3971  whichWay_ = saveWhichWay;
3972  xyBranchValue_ = saveXyBranchValue;
3973  chosen_ = saveChosen;
3974  return value;
3975}
3976OsiBiLinearBranchingObject::OsiBiLinearBranchingObject()
3977  :OsiTwoWayBranchingObject(),
3978   chosen_(0)
3979{
3980}
3981
3982// Useful constructor
3983OsiBiLinearBranchingObject::OsiBiLinearBranchingObject (OsiSolverInterface * solver,
3984                                                        const OsiBiLinear * set,
3985                                                        int way ,
3986                                                        double separator,
3987                                                        int chosen)
3988  :OsiTwoWayBranchingObject(solver,set,way,separator),
3989   chosen_(chosen)
3990{
3991  assert (chosen_>=0&&chosen_<2);
3992}
3993
3994// Copy constructor
3995OsiBiLinearBranchingObject::OsiBiLinearBranchingObject ( const OsiBiLinearBranchingObject & rhs) 
3996  :OsiTwoWayBranchingObject(rhs),
3997   chosen_(rhs.chosen_)
3998{
3999}
4000
4001// Assignment operator
4002OsiBiLinearBranchingObject & 
4003OsiBiLinearBranchingObject::operator=( const OsiBiLinearBranchingObject& rhs)
4004{
4005  if (this != &rhs) {
4006    OsiTwoWayBranchingObject::operator=(rhs);
4007    chosen_ = rhs.chosen_;
4008  }
4009  return *this;
4010}
4011OsiBranchingObject * 
4012OsiBiLinearBranchingObject::clone() const
4013{ 
4014  return (new OsiBiLinearBranchingObject(*this));
4015}
4016
4017
4018// Destructor
4019OsiBiLinearBranchingObject::~OsiBiLinearBranchingObject ()
4020{
4021}
4022double
4023OsiBiLinearBranchingObject::branch(OsiSolverInterface * solver)
4024{
4025  const OsiBiLinear * set =
4026    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
4027  assert (set);
4028  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
4029  branchIndex_++;
4030  set->newBounds(solver, way, chosen_, value_);
4031  return 0.0;
4032}
4033/* Return true if branch should only bound variables
4034 */
4035bool 
4036OsiBiLinearBranchingObject::boundBranch() const 
4037{ 
4038  const OsiBiLinear * set =
4039    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
4040  assert (set);
4041  return (set->branchingStrategy()&4)!=0;
4042}
4043// Print what would happen 
4044void
4045OsiBiLinearBranchingObject::print(const OsiSolverInterface * solver)
4046{
4047  const OsiBiLinear * set =
4048    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
4049  assert (set);
4050  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
4051  int iColumn = (chosen_==1) ? set->xColumn() : set->yColumn();
4052  printf("OsiBiLinear would branch %s on %c variable %d from value %g\n",
4053         (way<0) ? "down" : "up",
4054         (chosen_==0) ? 'X' : 'Y', iColumn, value_);
4055}
4056// Default Constructor
4057OsiBiLinearEquality::OsiBiLinearEquality ()
4058  : OsiBiLinear(),
4059    numberPoints_(0)
4060{
4061}
4062
4063// Useful constructor
4064OsiBiLinearEquality::OsiBiLinearEquality (OsiSolverInterface * solver, int xColumn,
4065                          int yColumn, int xyRow, double rhs,
4066                                          double xMesh)
4067  : OsiBiLinear(),
4068    numberPoints_(0)
4069{
4070  double xB[2];
4071  double yB[2];
4072  const double * lower = solver->getColLower();
4073  const double * upper = solver->getColUpper();
4074  xColumn_ = xColumn;
4075  yColumn_ = yColumn;
4076  xyRow_ = xyRow;
4077  coefficient_ = rhs;
4078  xB[0]=lower[xColumn_];
4079  xB[1]=upper[xColumn_];
4080  yB[0]=lower[yColumn_];
4081  yB[1]=upper[yColumn_];
4082  if (xB[1]*yB[1]<coefficient_+1.0e-12||
4083      xB[0]*yB[0]>coefficient_-1.0e-12) {
4084    printf("infeasible row - reformulate\n");
4085    abort();
4086  }
4087  // reduce range of x if possible
4088  if (yB[0]*xB[1]>coefficient_+1.0e12) {
4089    xB[1] = coefficient_/yB[0];
4090    solver->setColUpper(xColumn_,xB[1]);
4091  }
4092  if (yB[1]*xB[0]<coefficient_-1.0e12) {
4093    xB[0] = coefficient_/yB[1];
4094    solver->setColLower(xColumn_,xB[0]);
4095  }
4096  // See how many points
4097  numberPoints_ = (int) ((xB[1]-xB[0]+0.5*xMesh)/xMesh);
4098  // redo exactly
4099  xMeshSize_ = (xB[1]-xB[0])/((double) numberPoints_);
4100  numberPoints_++;
4101  //#define KEEPXY
4102#ifndef KEEPXY
4103  // Take out xyRow
4104  solver->setRowLower(xyRow_,0.0);
4105  solver->setRowUpper(xyRow_,0.0);
4106#else
4107  // make >=
4108  solver->setRowLower(xyRow_,coefficient_-0.05);
4109  solver->setRowUpper(xyRow_,COIN_DBL_MAX);
4110#endif
4111  double rowLower[3];
4112  double rowUpper[3];
4113#ifndef KEEPXY
4114  double * columnLower = new double [numberPoints_];
4115  double * columnUpper = new double [numberPoints_];
4116  double * objective = new double [numberPoints_];
4117  CoinBigIndex *starts = new CoinBigIndex[numberPoints_+1];
4118  int * index = new int[3*numberPoints_];
4119  double * element = new double [3*numberPoints_];
4120#else
4121  double * columnLower = new double [numberPoints_+2];
4122  double * columnUpper = new double [numberPoints_+2];
4123  double * objective = new double [numberPoints_+2];
4124  CoinBigIndex *starts = new CoinBigIndex[numberPoints_+3];
4125  int * index = new int[4*numberPoints_+2];
4126  double * element = new double [4*numberPoints_+2];
4127#endif
4128  int i;
4129  starts[0]=0;
4130  // rows
4131  int numberRows = solver->getNumRows();
4132  // convexity
4133  rowLower[0]=1.0;
4134  rowUpper[0]=1.0;
4135  convexity_ = numberRows;
4136  starts[1]=0;
4137  // x
4138  rowLower[1]=0.0;
4139  rowUpper[1]=0.0;
4140  index[0]=xColumn_;
4141  element[0]=-1.0;
4142  xRow_ = numberRows+1;
4143  starts[2]=1;
4144  rowLower[2]=0.0;
4145  rowUpper[2]=0.0;
4146  index[1]=yColumn;
4147  element[1]=-1.0;
4148  yRow_ = numberRows+2;
4149  starts[3]=2;
4150  solver->addRows(3,starts,index,element,rowLower,rowUpper);
4151  int n=0;
4152  firstLambda_ = solver->getNumCols();
4153  double x = xB[0];
4154  assert(xColumn_!=yColumn_);
4155  for (i=0;i<numberPoints_;i++) {
4156    double y = coefficient_/x;
4157    columnLower[i]=0.0;
4158    columnUpper[i]=2.0;
4159    objective[i]=0.0;
4160    double value;
4161#ifdef KEEPXY
4162    // xy
4163    value=coefficient_;
4164    element[n]=value;
4165    index[n++]=xyRow_;
4166#endif
4167    // convexity
4168    value=1.0;
4169    element[n]=value;
4170    index[n++]=0+numberRows;
4171    // x
4172    value=x;
4173    if (fabs(value)<1.0e-19)
4174      value = 1.0e-19;
4175    element[n]=value;
4176    index[n++]=1+numberRows;
4177    // y
4178    value=y;
4179    if (fabs(value)<1.0e-19)
4180      value = 1.0e-19;
4181    element[n]=value;
4182    index[n++]=2+numberRows;
4183    starts[i+1]=n;
4184    x += xMeshSize_;
4185  }
4186#ifdef KEEPXY
4187  // costed slacks
4188  columnLower[numberPoints_]=0.0;
4189  columnUpper[numberPoints_]=xMeshSize_;
4190  objective[numberPoints_]=1.0e3;;
4191  // convexity
4192  element[n]=1.0;
4193  index[n++]=0+numberRows;
4194  starts[numberPoints_+1]=n;
4195  columnLower[numberPoints_+1]=0.0;
4196  columnUpper[numberPoints_+1]=xMeshSize_;
4197  objective[numberPoints_+1]=1.0e3;;
4198  // convexity
4199  element[n]=-1.0;
4200  index[n++]=0+numberRows;
4201  starts[numberPoints_+2]=n;
4202  solver->addCols(numberPoints_+2,starts,index,element,columnLower,columnUpper,objective);
4203#else
4204  solver->addCols(numberPoints_,starts,index,element,columnLower,columnUpper,objective);
4205#endif
4206  delete [] columnLower;
4207  delete [] columnUpper;
4208  delete [] objective;
4209  delete [] starts;
4210  delete [] index;
4211  delete [] element;
4212}
4213
4214// Copy constructor
4215OsiBiLinearEquality::OsiBiLinearEquality ( const OsiBiLinearEquality & rhs)
4216  :OsiBiLinear(rhs),
4217   numberPoints_(rhs.numberPoints_)
4218{
4219}
4220
4221// Clone
4222OsiObject *
4223OsiBiLinearEquality::clone() const
4224{
4225  return new OsiBiLinearEquality(*this);
4226}
4227
4228// Assignment operator
4229OsiBiLinearEquality & 
4230OsiBiLinearEquality::operator=( const OsiBiLinearEquality& rhs)
4231{
4232  if (this!=&rhs) {
4233    OsiBiLinear::operator=(rhs);
4234    numberPoints_ = rhs.numberPoints_;
4235  }
4236  return *this;
4237}
4238
4239// Destructor
4240OsiBiLinearEquality::~OsiBiLinearEquality ()
4241{
4242}
4243// Possible improvement
4244double 
4245OsiBiLinearEquality::improvement(const OsiSolverInterface * solver) const
4246{
4247  const double * pi = solver->getRowPrice();
4248  int i;
4249  const double * solution = solver->getColSolution();
4250  printf(" for x %d y %d - pi %g %g\n",xColumn_,yColumn_,pi[xRow_],pi[yRow_]);
4251  for (i=0;i<numberPoints_;i++) {
4252    if (fabs(solution[i+firstLambda_])>1.0e-7)
4253      printf("(%d %g) ",i,solution[i+firstLambda_]);
4254  }
4255  printf("\n");
4256  return 0.0;
4257}
4258/* change grid
4259   if type 0 then use solution and make finer
4260   if 1 then back to original
4261*/
4262double 
4263OsiBiLinearEquality::newGrid(OsiSolverInterface * solver, int type) const
4264{
4265  CoinPackedMatrix * matrix = solver->getMutableMatrixByCol(); 
4266  if (!matrix) {
4267    printf("Unable to modify matrix\n");
4268    abort();
4269  }
4270  double * element = matrix->getMutableElements();
4271  const int * row = matrix->getIndices();
4272  const CoinBigIndex * columnStart = matrix->getVectorStarts();
4273  //const int * columnLength = matrix->getVectorLengths();
4274  // get original bounds
4275  double xB[2];
4276  double yB[2];
4277  const double * lower = solver->getColLower();
4278  const double * upper = solver->getColUpper();
4279  xB[0]=lower[xColumn_];
4280  xB[1]=upper[xColumn_];
4281  assert (fabs((xB[1]-xB[0])-xMeshSize_*(numberPoints_-1))<1.0e-7);
4282  yB[0]=lower[yColumn_];
4283  yB[1]=upper[yColumn_];
4284  double mesh=0.0;
4285  int i;
4286  if (type==0) {
4287    const double * solution = solver->getColSolution();
4288    int first=-1;
4289    int last=-1;
4290    double xValue=0.0;
4291    double step=0.0;
4292    for (i=0;i<numberPoints_;i++) {
4293      int iColumn = i+firstLambda_;
4294      if (fabs(solution[iColumn])>1.0e-7) {
4295        int k=columnStart[iColumn]+1;
4296        xValue += element[k]*solution[iColumn];
4297        if (first==-1) {
4298          first = i;
4299          step = -element[k];
4300        } else {
4301          step += element[k];
4302        }
4303        last =i;
4304      }
4305    }
4306    if (last>first+1) {
4307      printf("not adjacent - presuming small djs\n");
4308    }
4309    // new step size
4310    assert (numberPoints_>2);
4311    step = CoinMax((1.5*step)/((double) (numberPoints_-1)),0.5*step);
4312    xB[0] = CoinMax(xB[0],xValue-0.5*step);
4313    xB[1] = CoinMin(xB[1],xValue+0.5*step);
4314    // and now divide these
4315    mesh = (xB[1]-xB[0])/((double) (numberPoints_-1));
4316  } else {
4317    // back to original
4318    mesh = xMeshSize_;
4319  }
4320  double x = xB[0];
4321  for (i=0;i<numberPoints_;i++) {
4322    int iColumn = i+firstLambda_;
4323    double y = coefficient_/x;
4324    //assert (columnLength[iColumn]==3); - could have cuts
4325    int k=columnStart[iColumn];
4326#ifdef KEEPXY
4327    // xy
4328    assert (row[k]==xyRow_);
4329    k++;
4330#endif
4331    assert (row[k]==convexity_);
4332    k++;
4333    double value;
4334    // x
4335    value=x;
4336    assert (row[k]==xRow_);
4337    assert (fabs(value)>1.0e-10);
4338    element[k++]=value;
4339    // y
4340    value=y;
4341    assert (row[k]==yRow_);
4342    assert (fabs(value)>1.0e-10);
4343    element[k++]=value;
4344    x += mesh;
4345  }
4346  return mesh;
4347}
4348/** Default Constructor
4349
4350  Equivalent to an unspecified binary variable.
4351*/
4352OsiSimpleFixedInteger::OsiSimpleFixedInteger ()
4353  : OsiSimpleInteger()
4354{
4355}
4356
4357/** Useful constructor
4358
4359  Loads actual upper & lower bounds for the specified variable.
4360*/
4361OsiSimpleFixedInteger::OsiSimpleFixedInteger (const OsiSolverInterface * solver, int iColumn)
4362  : OsiSimpleInteger(solver,iColumn)
4363{
4364}
4365
4366 
4367// Useful constructor - passed solver index and original bounds
4368OsiSimpleFixedInteger::OsiSimpleFixedInteger ( int iColumn, double lower, double upper)
4369  : OsiSimpleInteger(iColumn,lower,upper)
4370{
4371}
4372
4373// Useful constructor - passed simple integer
4374OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleInteger &rhs)
4375  : OsiSimpleInteger(rhs)
4376{
4377}
4378
4379// Copy constructor
4380OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleFixedInteger & rhs)
4381  :OsiSimpleInteger(rhs)
4382
4383{
4384}
4385
4386// Clone
4387OsiObject *
4388OsiSimpleFixedInteger::clone() const
4389{
4390  return new OsiSimpleFixedInteger(*this);
4391}
4392
4393// Assignment operator
4394OsiSimpleFixedInteger & 
4395OsiSimpleFixedInteger::operator=( const OsiSimpleFixedInteger& rhs)
4396{
4397  if (this!=&rhs) {
4398    OsiSimpleInteger::operator=(rhs);
4399  }
4400  return *this;
4401}
4402
4403// Destructor
4404OsiSimpleFixedInteger::~OsiSimpleFixedInteger ()
4405{
4406}
4407// Infeasibility - large is 0.5
4408double 
4409OsiSimpleFixedInteger::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
4410{
4411  double value = info->solution_[columnNumber_];
4412  value = CoinMax(value, info->lower_[columnNumber_]);
4413  value = CoinMin(value, info->upper_[columnNumber_]);
4414  double nearest = floor(value+(1.0-0.5));
4415  if (nearest>value) { 
4416    whichWay=1;
4417  } else {
4418    whichWay=0;
4419  }
4420  infeasibility_ = fabs(value-nearest);
4421  bool satisfied=false;
4422  if (infeasibility_<=info->integerTolerance_) {
4423    otherInfeasibility_ = 1.0;
4424    satisfied=true;
4425    if (info->lower_[columnNumber_]!=info->upper_[columnNumber_])
4426      infeasibility_ = 1.0e-5;
4427    else
4428      infeasibility_ = 0.0;
4429  } else if (info->defaultDual_<0.0) {
4430    otherInfeasibility_ = 1.0-infeasibility_;
4431  } else {
4432    const double * pi = info->pi_;
4433    const double * activity = info->rowActivity_;
4434    const double * lower = info->rowLower_;
4435    const double * upper = info->rowUpper_;
4436    const double * element = info->elementByColumn_;
4437    const int * row = info->row_;
4438    const CoinBigIndex * columnStart = info->columnStart_;
4439    const int * columnLength = info->columnLength_;
4440    double direction = info->direction_;
4441    double downMovement = value - floor(value);
4442    double upMovement = 1.0-downMovement;
4443    double valueP = info->objective_[columnNumber_]*direction;
4444    CoinBigIndex start = columnStart[columnNumber_];
4445    CoinBigIndex end = start + columnLength[columnNumber_];
4446    double upEstimate = 0.0;
4447    double downEstimate = 0.0;
4448    if (valueP>0.0)
4449      upEstimate = valueP*upMovement;
4450    else
4451      downEstimate -= valueP*downMovement;
4452    double tolerance = info->primalTolerance_;
4453    for (CoinBigIndex j=start;j<end;j++) {
4454      int iRow = row[j];
4455      if (lower[iRow]<-1.0e20) 
4456        assert (pi[iRow]<=1.0e-3);
4457      if (upper[iRow]>1.0e20) 
4458        assert (pi[iRow]>=-1.0e-3);
4459      valueP = pi[iRow]*direction;
4460      double el2 = element[j];
4461      double value2 = valueP*el2;
4462      double u=0.0;
4463      double d=0.0;
4464      if (value2>0.0)
4465        u = value2;
4466      else
4467        d = -value2;
4468      // if up makes infeasible then make at least default
4469      double newUp = activity[iRow] + upMovement*el2;
4470      if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance)
4471        u = CoinMax(u,info->defaultDual_);
4472      upEstimate += u*upMovement;
4473      // if down makes infeasible then make at least default
4474      double newDown = activity[iRow] - downMovement*el2;
4475      if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance)
4476        d = CoinMax(d,info->defaultDual_);
4477      downEstimate += d*downMovement;
4478    }
4479    if (downEstimate>=upEstimate) {
4480      infeasibility_ = CoinMax(1.0e-12,upEstimate);
4481      otherInfeasibility_ = CoinMax(1.0e-12,downEstimate);
4482      whichWay = 1;
4483    } else {
4484      infeasibility_ = CoinMax(1.0e-12,downEstimate);
4485      otherInfeasibility_ = CoinMax(1.0e-12,upEstimate);
4486      whichWay = 0;
4487    }
4488  }
4489  if (preferredWay_>=0&&!satisfied)
4490    whichWay = preferredWay_;
4491  whichWay_=whichWay;
4492  return infeasibility_;
4493}
4494// Creates a branching object
4495OsiBranchingObject * 
4496OsiSimpleFixedInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const 
4497{
4498  double value = info->solution_[columnNumber_];
4499  value = CoinMax(value, info->lower_[columnNumber_]);
4500  value = CoinMin(value, info->upper_[columnNumber_]);
4501  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
4502  double nearest = floor(value+0.5);
4503  double integerTolerance = info->integerTolerance_;
4504  if (fabs(value-nearest)<integerTolerance) {
4505    // adjust value
4506    if (nearest!=info->upper_[columnNumber_])
4507      value = nearest+2.0*integerTolerance;
4508    else
4509      value = nearest-2.0*integerTolerance;
4510  }
4511  OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way,
4512                                             value);
4513  return branch;
4514}
4515
4516#include <cstdlib>
4517#include <cstdio>
4518#include <cmath>
4519#include <cfloat>
4520#include <cassert>
4521#include <iostream>
4522//#define CGL_DEBUG 2
4523#include "CoinHelperFunctions.hpp"
4524#include "CoinPackedVector.hpp"
4525#include "OsiRowCutDebugger.hpp"
4526#include "CoinWarmStartBasis.hpp"
4527//#include "CglTemporary.hpp"
4528#include "CoinFinite.hpp"
4529//-------------------------------------------------------------------
4530// Generate Stored cuts
4531//-------------------------------------------------------------------
4532void 
4533CglTemporary::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
4534                             const CglTreeInfo info) const
4535{
4536  // Get basic problem information
4537  const double * solution = si.getColSolution();
4538  int numberRowCuts = cuts_.sizeRowCuts();
4539  for (int i=0;i<numberRowCuts;i++) {
4540    const OsiRowCut * rowCutPointer = cuts_.rowCutPtr(i);
4541    double violation = rowCutPointer->violated(solution);
4542    if (violation>=requiredViolation_) 
4543      cs.insert(*rowCutPointer);
4544  }
4545  // delete
4546  cuts_ = OsiCuts();
4547}
4548
4549//-------------------------------------------------------------------
4550// Default Constructor
4551//-------------------------------------------------------------------
4552CglTemporary::CglTemporary ()
4553:
4554CglStored()
4555{
4556}
4557
4558//-------------------------------------------------------------------
4559// Copy constructor
4560//-------------------------------------------------------------------
4561CglTemporary::CglTemporary (const CglTemporary & source) :
4562  CglStored(source)
4563{ 
4564}
4565
4566//-------------------------------------------------------------------
4567// Clone
4568//-------------------------------------------------------------------
4569CglCutGenerator *
4570CglTemporary::clone() const
4571{
4572  return new CglTemporary(*this);
4573}
4574
4575//-------------------------------------------------------------------
4576// Destructor
4577//-------------------------------------------------------------------
4578CglTemporary::~CglTemporary ()
4579{
4580}
4581
4582//----------------------------------------------------------------
4583// Assignment operator
4584//-------------------------------------------------------------------
4585CglTemporary &
4586CglTemporary::operator=(const CglTemporary& rhs)
4587{
4588  if (this != &rhs) {
4589    CglStored::operator=(rhs);
4590  }
4591  return *this;
4592}
4593#endif
Note: See TracBrowser for help on using the repository browser.