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

Last change on this file since 493 was 493, checked in by forrest, 13 years ago

for quadratic integer etc

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