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

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

mods

File size: 115.6 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 (yRow_<0)
3003    yLambda = xLambda;
3004#if 0
3005  if (fabs(x-xLambda)>1.0e-4||
3006      fabs(y-yLambda)>1.0e-4)
3007    printf("feasibleregion x %d given %g lambda %g y %d given %g lambda %g\n",
3008           xColumn_,x,xLambda,
3009           yColumn_,y,yLambda);
3010#endif
3011#endif
3012  double infeasibility=0.0;
3013  double distance;
3014  double steps;
3015  double xNew=x;
3016  if (xMeshSize_) {
3017    distance = x-xB[0];
3018    if (x<0.5*(xB[0]+xB[1])) {
3019      distance = x-xB[0];
3020      steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
3021      xNew = xB[0]+steps*xMeshSize_;
3022      assert (xNew<=xB[1]+xSatisfied_);
3023    } else {
3024      distance = xB[1]-x;
3025      steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
3026      xNew = xB[1]-steps*xMeshSize_;
3027      assert (xNew>=xB[0]-xSatisfied_);
3028    }
3029    if (xMeshSize_<1.0&&fabs(xNew-x)<=xSatisfied_) {
3030      double lo = CoinMax(xB[0],x-0.5*xSatisfied_);
3031      double up = CoinMin(xB[1],x+0.5*xSatisfied_);
3032      solver->setColLower(xColumn_,lo);
3033      solver->setColUpper(xColumn_,up);
3034    } else {
3035      infeasibility +=  fabs(xNew-x);
3036      solver->setColLower(xColumn_,xNew);
3037      solver->setColUpper(xColumn_,xNew);
3038    }
3039  }
3040  double yNew=y;
3041  if (yMeshSize_) {
3042    distance = y-yB[0];
3043    if (y<0.5*(yB[0]+yB[1])) {
3044      distance = y-yB[0];
3045      steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
3046      yNew = yB[0]+steps*yMeshSize_;
3047      assert (yNew<=yB[1]+ySatisfied_);
3048    } else {
3049      distance = yB[1]-y;
3050      steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
3051      yNew = yB[1]-steps*yMeshSize_;
3052      assert (yNew>=yB[0]-ySatisfied_);
3053    }
3054    if (yMeshSize_<1.0&&fabs(yNew-y)<=ySatisfied_) {
3055      double lo = CoinMax(yB[0],y-0.5*ySatisfied_);
3056      double up = CoinMin(yB[1],y+0.5*ySatisfied_);
3057      solver->setColLower(yColumn_,lo);
3058      solver->setColUpper(yColumn_,up);
3059    } else {
3060      infeasibility +=  fabs(yNew-y);
3061      solver->setColLower(yColumn_,yNew);
3062      solver->setColUpper(yColumn_,yNew);
3063    }
3064  } 
3065  if (0) {
3066    // temp
3067      solver->setColLower(xColumn_,x);
3068      solver->setColUpper(xColumn_,x);
3069      solver->setColLower(yColumn_,y);
3070      solver->setColUpper(yColumn_,y);
3071  }
3072  if ((branchingStrategy_&4)) {
3073    // fake to make correct
3074    double lambda[4];
3075    computeLambdas(solver,lambda);
3076    for (int j=0;j<4;j++) {
3077      int iColumn = firstLambda_+j;
3078      double value = lambda[j];
3079      solver->setColLower(iColumn,value);
3080      solver->setColUpper(iColumn,value);
3081    }
3082  }
3083  double xyTrue = xNew*yNew;
3084  double xyLambda = 0.0;
3085  for (j=0;j<4;j++) {
3086    int iX = j>>1;
3087    int iY = j&1;
3088    xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j];
3089  }
3090  infeasibility += fabs(xyTrue-xyLambda);
3091  return infeasibility;
3092}
3093// Returns true value of single xyRow coefficient
3094double 
3095OsiBiLinear::xyCoefficient(const double * solution) const
3096{
3097  // If another object has finer mesh ignore this
3098  if ((branchingStrategy_&8)!=0)
3099    return 0.0;
3100  double x = solution[xColumn_];
3101  double y = solution[yColumn_];
3102  //printf("x (%d,%g) y (%d,%g) x*y*coefficient %g\n",
3103  // xColumn_,x,yColumn_,y,x*y*coefficient_);
3104  return x*y*coefficient_;
3105}
3106
3107// Redoes data when sequence numbers change
3108void 
3109OsiBiLinear::resetSequenceEtc(int numberColumns, const int * originalColumns)
3110{
3111  int i;
3112  for (i=0;i<numberColumns;i++) {
3113    if (originalColumns[i]==firstLambda_)
3114      break;
3115  }
3116  if (i<numberColumns) {
3117    firstLambda_ = i;
3118    for (int j=0;j<4;j++) {
3119      assert (originalColumns[j+i]-firstLambda_==j);
3120    }
3121  } else {
3122    printf("lost set\n");
3123    abort();
3124  }
3125  // rows will be out anyway
3126  abort();
3127}
3128
3129// Creates a branching object
3130OsiBranchingObject * 
3131OsiBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
3132{
3133  // create object
3134  OsiBranchingObject * branch;
3135  assert (chosen_==0||chosen_==1);
3136  //if (chosen_==0)
3137  //assert (xyBranchValue_>=info->lower_[xColumn_]&&xyBranchValue_<=info->upper_[xColumn_]);
3138  //else
3139  //assert (xyBranchValue_>=info->lower_[yColumn_]&&xyBranchValue_<=info->upper_[yColumn_]);
3140  branch = new OsiBiLinearBranchingObject(solver,this,way,xyBranchValue_,chosen_);
3141  return branch;
3142}
3143// Does work of branching
3144void 
3145OsiBiLinear::newBounds(OsiSolverInterface * solver, int way, short xOrY, double separator) const
3146{
3147  int iColumn;
3148  double mesh;
3149  double satisfied;
3150  if (xOrY==0) {
3151    iColumn=xColumn_;
3152    mesh=xMeshSize_;
3153    satisfied=xSatisfied_;
3154  } else {
3155    iColumn=yColumn_;
3156    mesh=yMeshSize_;
3157    satisfied=ySatisfied_;
3158  }
3159  const double * columnLower = solver->getColLower();
3160  const double * columnUpper = solver->getColUpper();
3161  double lower = columnLower[iColumn];
3162  double distance;
3163  double steps;
3164  double zNew=separator;
3165  distance = separator-lower;
3166  assert (mesh);
3167  steps = floor ((distance+0.5*mesh)/mesh);
3168  if (mesh<1.0)
3169    zNew = lower+steps*mesh;
3170  if (zNew>columnUpper[iColumn]-satisfied)
3171    zNew = 0.5*(columnUpper[iColumn]-lower);
3172  double oldUpper = columnUpper[iColumn] ;
3173  double oldLower = columnLower[iColumn] ;
3174#ifndef NDEBUG
3175  int nullChange=0;
3176#endif
3177  if (way<0) {
3178    if (zNew>separator&&mesh<1.0)
3179      zNew -= mesh;
3180    double oldUpper = columnUpper[iColumn] ;
3181    if (zNew+satisfied>=oldUpper)
3182      zNew =0.5*(oldUpper+oldLower);
3183    if (mesh==1.0)
3184      zNew = floor(separator);
3185#ifndef NDEBUG
3186    if (oldUpper<zNew+1.0e-8)
3187      nullChange=-1;
3188#endif
3189    solver->setColUpper(iColumn,zNew);
3190  } else {
3191    if (zNew<separator&&mesh<1.0)
3192      zNew += mesh;
3193    if (zNew-satisfied<=oldLower)
3194      zNew =0.5*(oldUpper+oldLower);
3195    if (mesh==1.0)
3196      zNew = ceil(separator);
3197#ifndef NDEBUG
3198    if (oldLower>zNew-1.0e-8)
3199      nullChange=1;
3200#endif
3201    solver->setColLower(iColumn,zNew);
3202  }
3203  if ((branchingStrategy_&4)!=0
3204      &&columnLower[xColumn_]==columnUpper[xColumn_]
3205      &&columnLower[yColumn_]==columnUpper[yColumn_]) {
3206    // fake to make correct
3207    double lambda[4];
3208    computeLambdas(solver,lambda);
3209    for (int j=0;j<4;j++) {
3210      int iColumn = firstLambda_+j;
3211      double value = lambda[j];
3212#ifndef NDEBUG
3213      if (fabs(value-columnLower[iColumn])>1.0e-5||
3214          fabs(value-columnUpper[iColumn])>1.0e-5)
3215        nullChange=0;
3216#endif
3217      solver->setColLower(iColumn,value);
3218      solver->setColUpper(iColumn,value);
3219    }
3220  }
3221#ifndef NDEBUG
3222  if (nullChange)
3223    printf("null change on column%s %d - bounds %g,%g\n",nullChange>0 ? "Lower" : "Upper",
3224           iColumn,oldLower,oldUpper);
3225#endif
3226#if 0
3227  // always free up lambda
3228  for (int i=firstLambda_;i<firstLambda_+4;i++) {
3229    solver->setColLower(i,0.0);
3230    solver->setColUpper(i,2.0);
3231  }
3232#endif
3233  double xB[3];
3234  xB[0]=columnLower[xColumn_];
3235  xB[1]=columnUpper[xColumn_];
3236  double yB[3];
3237  yB[0]=columnLower[yColumn_];
3238  yB[1]=columnUpper[yColumn_];
3239  if (false&&(branchingStrategy_&4)!=0&&yRow_>=0&&
3240      xMeshSize_==1.0&&yMeshSize_==1.0) {
3241    if ((xB[1]-xB[0])*(yB[1]-yB[0])<40) {
3242      // try looking at all solutions
3243      double lower[4];
3244      double upper[4];
3245      double lambda[4];
3246      int i;
3247      double lowerLambda[4];
3248      double upperLambda[4];
3249      for (i=0;i<4;i++) {
3250        lower[i] = CoinMax(0.0,columnLower[firstLambda_+i]);
3251        upper[i] = CoinMin(1.0,columnUpper[firstLambda_+i]);
3252        lowerLambda[i]=1.0;
3253        upperLambda[i]=0.0;
3254      }
3255      // get coefficients
3256      double xybar[4];
3257      getCoefficients(solver,xB,yB,xybar);
3258      double x,y;
3259      for (x=xB[0];x<=xB[1];x++) {
3260        xB[2]=x;
3261        for (y=yB[0];y<=yB[1];y++) {
3262          yB[2]=y;
3263          computeLambdas(xB, yB, xybar,lambda);
3264          for (i=0;i<4;i++) {
3265            lowerLambda[i] = CoinMin(lowerLambda[i],lambda[i]);
3266            upperLambda[i] = CoinMax(upperLambda[i],lambda[i]);
3267          }
3268        }
3269      }
3270      double change=0.0;;
3271      for (i=0;i<4;i++) {
3272        if (lowerLambda[i]>lower[i]+1.0e-12) {
3273          solver->setColLower(firstLambda_+i,lowerLambda[i]);
3274          change += lowerLambda[i]-lower[i];
3275        }
3276        if (upperLambda[i]<upper[i]-1.0e-12) {
3277          solver->setColUpper(firstLambda_+i,upperLambda[i]);
3278          change -= upperLambda[i]-upper[i];
3279        }
3280      }
3281      if (change>1.0e-5)
3282        printf("change of %g\n",change);
3283    }
3284  }
3285  if (boundType_) {
3286    assert (!xMeshSize_||!yMeshSize_);
3287    if (xMeshSize_) {
3288      // can tighten bounds on y
3289      if ((boundType_&1)!=0) {
3290        if (xB[0]*yB[1]>coefficient_) {
3291          // tighten upper bound on y
3292          solver->setColUpper(yColumn_,coefficient_/xB[0]);
3293        }
3294      }
3295      if ((boundType_&2)!=0) {
3296        if (xB[1]*yB[0]<coefficient_) {
3297          // tighten lower bound on y
3298          solver->setColLower(yColumn_,coefficient_/xB[1]);
3299        }
3300      }
3301    } else {
3302      // can tighten bounds on x
3303      if ((boundType_&1)!=0) {
3304        if (yB[0]*xB[1]>coefficient_) {
3305          // tighten upper bound on x
3306          solver->setColUpper(xColumn_,coefficient_/yB[0]);
3307        }
3308      }
3309      if ((boundType_&2)!=0) {
3310        if (yB[1]*xB[0]<coefficient_) {
3311          // tighten lower bound on x
3312          solver->setColLower(xColumn_,coefficient_/yB[1]);
3313        }
3314      }
3315    }
3316  }
3317}
3318// Compute lambdas if coefficients not changing
3319void 
3320OsiBiLinear::computeLambdas(const OsiSolverInterface * solver,double lambda[4]) const
3321{
3322  // fix so correct
3323  double xB[3],yB[3];
3324  double xybar[4];
3325  getCoefficients(solver,xB,yB,xybar);
3326  double x,y;
3327  x=solver->getColLower()[xColumn_];
3328  assert(x==solver->getColUpper()[xColumn_]);
3329  xB[2]=x;
3330  y=solver->getColLower()[yColumn_];
3331  assert(y==solver->getColUpper()[yColumn_]);
3332  yB[2]=y;
3333  computeLambdas(xB, yB, xybar,lambda);
3334  assert (xyRow_>=0);
3335}
3336// Get LU coefficients from matrix
3337void 
3338OsiBiLinear::getCoefficients(const OsiSolverInterface * solver,double xB[2], double yB[2],
3339                             double xybar[4]) const
3340{
3341  const CoinPackedMatrix * matrix = solver->getMatrixByCol();
3342  const double * element = matrix->getElements();
3343  const double * objective = solver->getObjCoefficients();
3344  const int * row = matrix->getIndices();
3345  const CoinBigIndex * columnStart = matrix->getVectorStarts();
3346  const int * columnLength = matrix->getVectorLengths();
3347  // order is LxLy, LxUy, UxLy and UxUy
3348  int j;
3349  double multiplier = (boundType_==0) ? 1.0/coefficient_ : 1.0;
3350  for (j=0;j<4;j++) {
3351    int iColumn = firstLambda_+j;
3352    int iStart = columnStart[iColumn];
3353    int iEnd = iStart + columnLength[iColumn];
3354    int k=iStart;
3355    double x=0.0;
3356    double y=0.0;
3357    xybar[j]=0.0;
3358    for (;k<iEnd;k++) {
3359      if (xRow_==row[k])
3360        x = element[k];
3361      if (yRow_==row[k])
3362        y = element[k];
3363      if (xyRow_==row[k])
3364        xybar[j] = element[k]*multiplier;
3365    }
3366    if (yRow_<0)
3367      y=x;
3368    if (xyRow_<0)
3369      xybar[j] = objective[iColumn]*multiplier;
3370    assert (fabs(xybar[j]-x*y)<1.0e-4);
3371    if (j==0)
3372      xB[0]=x;
3373    else if (j==1)
3374      yB[1]=y;
3375    else if (j==2)
3376      yB[0]=y;
3377    else if (j==3)
3378      xB[1]=x;
3379  }
3380}
3381// Compute lambdas (third entry in each .B is current value)
3382double
3383OsiBiLinear::computeLambdas(const double xB[3], const double yB[3], const double xybar[4], double lambda[4]) const
3384{
3385  // fake to make correct
3386  double x = xB[2];
3387  double y = yB[2];
3388  // order is LxLy, LxUy, UxLy and UxUy
3389  // l0 + l1 = this
3390  double rhs1 = (xB[1]-x)/(xB[1]-xB[0]);
3391  // l0 + l2 = this
3392  double rhs2 = (yB[1]-y)/(yB[1]-yB[0]);
3393  // For xy (taking out l3)
3394  double rhs3 = xB[1]*yB[1]-x*y;
3395  double a0 = xB[1]*yB[1] - xB[0]*yB[0];
3396  double a1 = xB[1]*yB[1] - xB[0]*yB[1];
3397  double a2 = xB[1]*yB[1] - xB[1]*yB[0];
3398  // divide through to get l0 coefficient
3399  rhs3 /= a0;
3400  a1 /= a0;
3401  a2 /= a0;
3402  // subtract out l0
3403  double b[2][2];
3404  double rhs[2];
3405  // first for l1 and l2
3406  b[0][0] = 1.0-a1;
3407  b[0][1] = -a2;
3408  rhs[0] = rhs1 - rhs3;
3409  // second
3410  b[1][0] = -a1;
3411  b[1][1] = 1.0-a2;
3412  rhs[1] = rhs2 - rhs3;
3413  if (fabs(b[0][0])>fabs(b[0][1])) {
3414    double sub = b[1][0]/b[0][0];
3415    b[1][1] -= sub*b[0][1];
3416    rhs[1] -= sub*rhs[0];
3417    assert (fabs(b[1][1])>1.0e-12);
3418    lambda[2] = rhs[1]/b[1][1];
3419    lambda[0] = rhs2 - lambda[2];
3420    lambda[1] = rhs1 - lambda[0];
3421  } else {
3422    double sub = b[1][1]/b[0][1];
3423    b[1][0] -= sub*b[0][0];
3424    rhs[1] -= sub*rhs[0];
3425    assert (fabs(b[1][0])>1.0e-12);
3426    lambda[1] = rhs[1]/b[1][0];
3427    lambda[0] = rhs1 - lambda[1];
3428    lambda[2] = rhs2 - lambda[0];
3429  }
3430  lambda[3] = 1.0- (lambda[0]+lambda[1]+lambda[2]);
3431  double infeasibility=0.0;
3432  double xy=0.0;
3433  for (int j=0;j<4;j++) {
3434    double value = lambda[j];
3435    if (value>1.0) {
3436      infeasibility += value-1.0;
3437      value=1.0;
3438    }
3439    if (value<0.0) {
3440      infeasibility -= value;
3441      value=0.0;
3442    }
3443    lambda[j]=value;
3444    xy += xybar[j]*value;
3445  }
3446  assert (fabs(xy-x*y)<1.0e-4);
3447  return infeasibility;
3448}
3449// Updates coefficients
3450int
3451OsiBiLinear::updateCoefficients(const double * lower, const double * upper, double * objective, 
3452                                CoinPackedMatrix * matrix, CoinWarmStartBasis * basis) const
3453{
3454  // Return if no updates
3455  if ((branchingStrategy_&4)!=0)
3456    return 0;
3457  int numberUpdated=0;
3458  double * element = matrix->getMutableElements();
3459  const int * row = matrix->getIndices();
3460  const CoinBigIndex * columnStart = matrix->getVectorStarts();
3461  //const int * columnLength = matrix->getVectorLengths();
3462  // order is LxLy, LxUy, UxLy and UxUy
3463  double xB[2];
3464  double yB[2];
3465  xB[0]=lower[xColumn_];
3466  xB[1]=upper[xColumn_];
3467  yB[0]=lower[yColumn_];
3468  yB[1]=upper[yColumn_];
3469  //printf("x %d (%g,%g) y %d (%g,%g)\n",
3470  // xColumn_,xB[0],xB[1],
3471  // yColumn_,yB[0],yB[1]);
3472  CoinWarmStartBasis::Status status[4];
3473  int numStruct = basis ? basis->getNumStructural()-firstLambda_ : 0;
3474  double coefficient = (boundType_==0) ? coefficient_ : 1.0;
3475  for (int j=0;j<4;j++) {
3476    status[j]=(j<numStruct) ? basis->getStructStatus(j+firstLambda_) : CoinWarmStartBasis::atLowerBound;
3477    int iX = j>>1;
3478    double x = xB[iX];
3479    int iY = j&1;
3480    double y = yB[iY];
3481    CoinBigIndex k = columnStart[j+firstLambda_];
3482    double value;
3483    // xy
3484    value=coefficient*x*y;
3485    if (xyRow_>=0) {
3486      assert (row[k]==xyRow_);
3487#if BI_PRINT > 1
3488      printf("j %d xy (%d,%d) coeff from %g to %g\n",j,xColumn_,yColumn_,element[k],value);
3489#endif
3490      element[k++]=value;
3491    } else {
3492      // objective
3493      objective[j+firstLambda_]=value;
3494    }
3495    numberUpdated++;
3496    // convexity
3497    assert (row[k]==convexity_);
3498    k++;
3499    // x
3500    value=x;
3501#if BI_PRINT > 1
3502    printf("j %d x (%d) coeff from %g to %g\n",j,xColumn_,element[k],value);
3503#endif
3504    assert (row[k]==xRow_);
3505    element[k++]=value;
3506    numberUpdated++;
3507    if (yRow_>=0) {
3508      // y
3509      value=y;
3510#if BI_PRINT > 1
3511      printf("j %d y (%d) coeff from %g to %g\n",j,yColumn_,element[k],value);
3512#endif
3513      assert (row[k]==yRow_);
3514      element[k++]=value;
3515      numberUpdated++;
3516    }
3517  }
3518 
3519  if (xB[0]==xB[1]) {
3520    if (yB[0]==yB[1]) {
3521      // only one basic
3522      bool first=true;
3523      for (int j=0;j<4;j++) {
3524        if (status[j]==CoinWarmStartBasis::basic) {
3525          if (first) {
3526            first=false;
3527          } else {
3528            basis->setStructStatus(j+firstLambda_,CoinWarmStartBasis::atLowerBound);
3529#if BI_PRINT
3530            printf("zapping %d (x=%d,y=%d)\n",j,xColumn_,yColumn_);
3531#endif
3532          }
3533        }
3534      }
3535    } else {
3536      if (status[0]==CoinWarmStartBasis::basic&&
3537          status[2]==CoinWarmStartBasis::basic) {
3538        basis->setStructStatus(2+firstLambda_,CoinWarmStartBasis::atLowerBound);
3539#if BI_PRINT
3540        printf("zapping %d (x=%d,y=%d)\n",2,xColumn_,yColumn_);
3541#endif
3542      }
3543      if (status[1]==CoinWarmStartBasis::basic&&
3544          status[3]==CoinWarmStartBasis::basic) {
3545        basis->setStructStatus(3+firstLambda_,CoinWarmStartBasis::atLowerBound);
3546#if BI_PRINT
3547        printf("zapping %d (x=%d,y=%d)\n",3,xColumn_,yColumn_);
3548#endif
3549      }
3550    }
3551  } else if (yB[0]==yB[1]) {
3552    if (status[0]==CoinWarmStartBasis::basic&&
3553        status[1]==CoinWarmStartBasis::basic) {
3554      basis->setStructStatus(1+firstLambda_,CoinWarmStartBasis::atLowerBound);
3555#if BI_PRINT
3556      printf("zapping %d (x=%d,y=%d)\n",1,xColumn_,yColumn_);
3557#endif
3558    }
3559    if (status[2]==CoinWarmStartBasis::basic&&
3560        status[3]==CoinWarmStartBasis::basic) {
3561      basis->setStructStatus(3+firstLambda_,CoinWarmStartBasis::atLowerBound);
3562#if BI_PRINT
3563      printf("zapping %d (x=%d,y=%d)\n",3,xColumn_,yColumn_);
3564#endif
3565    }
3566  }
3567  return numberUpdated;
3568}
3569// This does NOT set mutable stuff
3570double 
3571OsiBiLinear::checkInfeasibility(const OsiBranchingInformation * info) const
3572{
3573  // If another object has finer mesh ignore this
3574  if ((branchingStrategy_&8)!=0)
3575    return 0.0;
3576  int way;
3577  double saveInfeasibility = infeasibility_;
3578  int saveWhichWay = whichWay_;
3579  double saveXyBranchValue = xyBranchValue_;
3580  short saveChosen = chosen_;
3581  double value = infeasibility(info,way);
3582  infeasibility_ = saveInfeasibility;
3583  whichWay_ = saveWhichWay;
3584  xyBranchValue_ = saveXyBranchValue;
3585  chosen_ = saveChosen;
3586  return value;
3587}
3588OsiBiLinearBranchingObject::OsiBiLinearBranchingObject()
3589  :OsiTwoWayBranchingObject(),
3590   chosen_(0)
3591{
3592}
3593
3594// Useful constructor
3595OsiBiLinearBranchingObject::OsiBiLinearBranchingObject (OsiSolverInterface * solver,
3596                                                        const OsiBiLinear * set,
3597                                                        int way ,
3598                                                        double separator,
3599                                                        int chosen)
3600  :OsiTwoWayBranchingObject(solver,set,way,separator),
3601   chosen_(chosen)
3602{
3603  assert (chosen_>=0&&chosen_<2);
3604}
3605
3606// Copy constructor
3607OsiBiLinearBranchingObject::OsiBiLinearBranchingObject ( const OsiBiLinearBranchingObject & rhs) 
3608  :OsiTwoWayBranchingObject(rhs),
3609   chosen_(rhs.chosen_)
3610{
3611}
3612
3613// Assignment operator
3614OsiBiLinearBranchingObject & 
3615OsiBiLinearBranchingObject::operator=( const OsiBiLinearBranchingObject& rhs)
3616{
3617  if (this != &rhs) {
3618    OsiTwoWayBranchingObject::operator=(rhs);
3619    chosen_ = rhs.chosen_;
3620  }
3621  return *this;
3622}
3623OsiBranchingObject * 
3624OsiBiLinearBranchingObject::clone() const
3625{ 
3626  return (new OsiBiLinearBranchingObject(*this));
3627}
3628
3629
3630// Destructor
3631OsiBiLinearBranchingObject::~OsiBiLinearBranchingObject ()
3632{
3633}
3634double
3635OsiBiLinearBranchingObject::branch(OsiSolverInterface * solver)
3636{
3637  const OsiBiLinear * set =
3638    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
3639  assert (set);
3640  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
3641  branchIndex_++;
3642  set->newBounds(solver, way, chosen_, value_);
3643  return 0.0;
3644}
3645/* Return true if branch should only bound variables
3646 */
3647bool 
3648OsiBiLinearBranchingObject::boundBranch() const 
3649{ 
3650  const OsiBiLinear * set =
3651    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
3652  assert (set);
3653  return (set->branchingStrategy()&4)!=0;
3654}
3655// Print what would happen 
3656void
3657OsiBiLinearBranchingObject::print(const OsiSolverInterface * solver)
3658{
3659  const OsiBiLinear * set =
3660    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
3661  assert (set);
3662  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
3663  int iColumn = (chosen_==1) ? set->xColumn() : set->yColumn();
3664  printf("OsiBiLinear would branch %s on %c variable %d from value %g\n",
3665         (way<0) ? "down" : "up",
3666         (chosen_==0) ? 'X' : 'Y', iColumn, value_);
3667}
3668// Default Constructor
3669OsiBiLinearEquality::OsiBiLinearEquality ()
3670  : OsiBiLinear(),
3671    numberPoints_(0)
3672{
3673}
3674
3675// Useful constructor
3676OsiBiLinearEquality::OsiBiLinearEquality (OsiSolverInterface * solver, int xColumn,
3677                          int yColumn, int xyRow, double rhs,
3678                                          double xMesh)
3679  : OsiBiLinear(),
3680    numberPoints_(0)
3681{
3682  double xB[2];
3683  double yB[2];
3684  const double * lower = solver->getColLower();
3685  const double * upper = solver->getColUpper();
3686  xColumn_ = xColumn;
3687  yColumn_ = yColumn;
3688  xyRow_ = xyRow;
3689  coefficient_ = rhs;
3690  xB[0]=lower[xColumn_];
3691  xB[1]=upper[xColumn_];
3692  yB[0]=lower[yColumn_];
3693  yB[1]=upper[yColumn_];
3694  if (xB[1]*yB[1]<coefficient_+1.0e-12||
3695      xB[0]*yB[0]>coefficient_-1.0e-12) {
3696    printf("infeasible row - reformulate\n");
3697    abort();
3698  }
3699  // reduce range of x if possible
3700  if (yB[0]*xB[1]>coefficient_+1.0e12) {
3701    xB[1] = coefficient_/yB[0];
3702    solver->setColUpper(xColumn_,xB[1]);
3703  }
3704  if (yB[1]*xB[0]<coefficient_-1.0e12) {
3705    xB[0] = coefficient_/yB[1];
3706    solver->setColLower(xColumn_,xB[0]);
3707  }
3708  // See how many points
3709  numberPoints_ = (int) ((xB[1]-xB[0]+0.5*xMesh)/xMesh);
3710  // redo exactly
3711  xMeshSize_ = (xB[1]-xB[0])/((double) numberPoints_);
3712  numberPoints_++;
3713  //#define KEEPXY
3714#ifndef KEEPXY
3715  // Take out xyRow
3716  solver->setRowLower(xyRow_,0.0);
3717  solver->setRowUpper(xyRow_,0.0);
3718#else
3719  // make >=
3720  solver->setRowLower(xyRow_,coefficient_-0.05);
3721  solver->setRowUpper(xyRow_,COIN_DBL_MAX);
3722#endif
3723  double rowLower[3];
3724  double rowUpper[3];
3725#ifndef KEEPXY
3726  double * columnLower = new double [numberPoints_];
3727  double * columnUpper = new double [numberPoints_];
3728  double * objective = new double [numberPoints_];
3729  CoinBigIndex *starts = new CoinBigIndex[numberPoints_+1];
3730  int * index = new int[3*numberPoints_];
3731  double * element = new double [3*numberPoints_];
3732#else
3733  double * columnLower = new double [numberPoints_+2];
3734  double * columnUpper = new double [numberPoints_+2];
3735  double * objective = new double [numberPoints_+2];
3736  CoinBigIndex *starts = new CoinBigIndex[numberPoints_+3];
3737  int * index = new int[4*numberPoints_+2];
3738  double * element = new double [4*numberPoints_+2];
3739#endif
3740  int i;
3741  starts[0]=0;
3742  // rows
3743  int numberRows = solver->getNumRows();
3744  // convexity
3745  rowLower[0]=1.0;
3746  rowUpper[0]=1.0;
3747  convexity_ = numberRows;
3748  starts[1]=0;
3749  // x
3750  rowLower[1]=0.0;
3751  rowUpper[1]=0.0;
3752  index[0]=xColumn_;
3753  element[0]=-1.0;
3754  xRow_ = numberRows+1;
3755  starts[2]=1;
3756  rowLower[2]=0.0;
3757  rowUpper[2]=0.0;
3758  index[1]=yColumn;
3759  element[1]=-1.0;
3760  yRow_ = numberRows+2;
3761  starts[3]=2;
3762  solver->addRows(3,starts,index,element,rowLower,rowUpper);
3763  int n=0;
3764  firstLambda_ = solver->getNumCols();
3765  double x = xB[0];
3766  assert(xColumn_!=yColumn_);
3767  for (i=0;i<numberPoints_;i++) {
3768    double y = coefficient_/x;
3769    columnLower[i]=0.0;
3770    columnUpper[i]=2.0;
3771    objective[i]=0.0;
3772    double value;
3773#ifdef KEEPXY
3774    // xy
3775    value=coefficient_;
3776    element[n]=value;
3777    index[n++]=xyRow_;
3778#endif
3779    // convexity
3780    value=1.0;
3781    element[n]=value;
3782    index[n++]=0+numberRows;
3783    // x
3784    value=x;
3785    if (fabs(value)<1.0e-19)
3786      value = 1.0e-19;
3787    element[n]=value;
3788    index[n++]=1+numberRows;
3789    // y
3790    value=y;
3791    if (fabs(value)<1.0e-19)
3792      value = 1.0e-19;
3793    element[n]=value;
3794    index[n++]=2+numberRows;
3795    starts[i+1]=n;
3796    x += xMeshSize_;
3797  }
3798#ifdef KEEPXY
3799  // costed slacks
3800  columnLower[numberPoints_]=0.0;
3801  columnUpper[numberPoints_]=xMeshSize_;
3802  objective[numberPoints_]=1.0e3;;
3803  // convexity
3804  element[n]=1.0;
3805  index[n++]=0+numberRows;
3806  starts[numberPoints_+1]=n;
3807  columnLower[numberPoints_+1]=0.0;
3808  columnUpper[numberPoints_+1]=xMeshSize_;
3809  objective[numberPoints_+1]=1.0e3;;
3810  // convexity
3811  element[n]=-1.0;
3812  index[n++]=0+numberRows;
3813  starts[numberPoints_+2]=n;
3814  solver->addCols(numberPoints_+2,starts,index,element,columnLower,columnUpper,objective);
3815#else
3816  solver->addCols(numberPoints_,starts,index,element,columnLower,columnUpper,objective);
3817#endif
3818  delete [] columnLower;
3819  delete [] columnUpper;
3820  delete [] objective;
3821  delete [] starts;
3822  delete [] index;
3823  delete [] element;
3824}
3825
3826// Copy constructor
3827OsiBiLinearEquality::OsiBiLinearEquality ( const OsiBiLinearEquality & rhs)
3828  :OsiBiLinear(rhs),
3829   numberPoints_(rhs.numberPoints_)
3830{
3831}
3832
3833// Clone
3834OsiObject *
3835OsiBiLinearEquality::clone() const
3836{
3837  return new OsiBiLinearEquality(*this);
3838}
3839
3840// Assignment operator
3841OsiBiLinearEquality & 
3842OsiBiLinearEquality::operator=( const OsiBiLinearEquality& rhs)
3843{
3844  if (this!=&rhs) {
3845    OsiBiLinear::operator=(rhs);
3846    numberPoints_ = rhs.numberPoints_;
3847  }
3848  return *this;
3849}
3850
3851// Destructor
3852OsiBiLinearEquality::~OsiBiLinearEquality ()
3853{
3854}
3855// Possible improvement
3856double 
3857OsiBiLinearEquality::improvement(const OsiSolverInterface * solver) const
3858{
3859  const double * pi = solver->getRowPrice();
3860  int i;
3861  const double * solution = solver->getColSolution();
3862  printf(" for x %d y %d - pi %g %g\n",xColumn_,yColumn_,pi[xRow_],pi[yRow_]);
3863  for (i=0;i<numberPoints_;i++) {
3864    if (fabs(solution[i+firstLambda_])>1.0e-7)
3865      printf("(%d %g) ",i,solution[i+firstLambda_]);
3866  }
3867  printf("\n");
3868  return 0.0;
3869}
3870/* change grid
3871   if type 0 then use solution and make finer
3872   if 1 then back to original
3873*/
3874double 
3875OsiBiLinearEquality::newGrid(OsiSolverInterface * solver, int type) const
3876{
3877  CoinPackedMatrix * matrix = solver->getMutableMatrixByCol(); 
3878  if (!matrix) {
3879    printf("Unable to modify matrix\n");
3880    abort();
3881  }
3882  double * element = matrix->getMutableElements();
3883  const int * row = matrix->getIndices();
3884  const CoinBigIndex * columnStart = matrix->getVectorStarts();
3885  //const int * columnLength = matrix->getVectorLengths();
3886  // get original bounds
3887  double xB[2];
3888  double yB[2];
3889  const double * lower = solver->getColLower();
3890  const double * upper = solver->getColUpper();
3891  xB[0]=lower[xColumn_];
3892  xB[1]=upper[xColumn_];
3893  assert (fabs((xB[1]-xB[0])-xMeshSize_*(numberPoints_-1))<1.0e-7);
3894  yB[0]=lower[yColumn_];
3895  yB[1]=upper[yColumn_];
3896  double mesh=0.0;
3897  int i;
3898  if (type==0) {
3899    const double * solution = solver->getColSolution();
3900    int first=-1;
3901    int last=-1;
3902    double xValue=0.0;
3903    double step=0.0;
3904    for (i=0;i<numberPoints_;i++) {
3905      int iColumn = i+firstLambda_;
3906      if (fabs(solution[iColumn])>1.0e-7) {
3907        int k=columnStart[iColumn]+1;
3908        xValue += element[k]*solution[iColumn];
3909        if (first==-1) {
3910          first = i;
3911          step = -element[k];
3912        } else {
3913          step += element[k];
3914        }
3915        last =i;
3916      }
3917    }
3918    if (last>first+1) {
3919      printf("not adjacent - presuming small djs\n");
3920    }
3921    // new step size
3922    assert (numberPoints_>2);
3923    step = CoinMax((1.5*step)/((double) (numberPoints_-1)),0.5*step);
3924    xB[0] = CoinMax(xB[0],xValue-0.5*step);
3925    xB[1] = CoinMin(xB[1],xValue+0.5*step);
3926    // and now divide these
3927    mesh = (xB[1]-xB[0])/((double) (numberPoints_-1));
3928  } else {
3929    // back to original
3930    mesh = xMeshSize_;
3931  }
3932  double x = xB[0];
3933  for (i=0;i<numberPoints_;i++) {
3934    int iColumn = i+firstLambda_;
3935    double y = coefficient_/x;
3936    //assert (columnLength[iColumn]==3); - could have cuts
3937    int k=columnStart[iColumn];
3938#ifdef KEEPXY
3939    // xy
3940    assert (row[k]==xyRow_);
3941    k++;
3942#endif
3943    assert (row[k]==convexity_);
3944    k++;
3945    double value;
3946    // x
3947    value=x;
3948    assert (row[k]==xRow_);
3949    assert (fabs(value)>1.0e-10);
3950    element[k++]=value;
3951    // y
3952    value=y;
3953    assert (row[k]==yRow_);
3954    assert (fabs(value)>1.0e-10);
3955    element[k++]=value;
3956    x += mesh;
3957  }
3958  return mesh;
3959}
3960/** Default Constructor
3961
3962  Equivalent to an unspecified binary variable.
3963*/
3964OsiSimpleFixedInteger::OsiSimpleFixedInteger ()
3965  : OsiSimpleInteger()
3966{
3967}
3968
3969/** Useful constructor
3970
3971  Loads actual upper & lower bounds for the specified variable.
3972*/
3973OsiSimpleFixedInteger::OsiSimpleFixedInteger (const OsiSolverInterface * solver, int iColumn)
3974  : OsiSimpleInteger(solver,iColumn)
3975{
3976}
3977
3978 
3979// Useful constructor - passed solver index and original bounds
3980OsiSimpleFixedInteger::OsiSimpleFixedInteger ( int iColumn, double lower, double upper)
3981  : OsiSimpleInteger(iColumn,lower,upper)
3982{
3983}
3984
3985// Useful constructor - passed simple integer
3986OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleInteger &rhs)
3987  : OsiSimpleInteger(rhs)
3988{
3989}
3990
3991// Copy constructor
3992OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleFixedInteger & rhs)
3993  :OsiSimpleInteger(rhs)
3994
3995{
3996}
3997
3998// Clone
3999OsiObject *
4000OsiSimpleFixedInteger::clone() const
4001{
4002  return new OsiSimpleFixedInteger(*this);
4003}
4004
4005// Assignment operator
4006OsiSimpleFixedInteger & 
4007OsiSimpleFixedInteger::operator=( const OsiSimpleFixedInteger& rhs)
4008{
4009  if (this!=&rhs) {
4010    OsiSimpleInteger::operator=(rhs);
4011  }
4012  return *this;
4013}
4014
4015// Destructor
4016OsiSimpleFixedInteger::~OsiSimpleFixedInteger ()
4017{
4018}
4019// Infeasibility - large is 0.5
4020double 
4021OsiSimpleFixedInteger::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
4022{
4023  double value = info->solution_[columnNumber_];
4024  value = CoinMax(value, info->lower_[columnNumber_]);
4025  value = CoinMin(value, info->upper_[columnNumber_]);
4026  double nearest = floor(value+(1.0-0.5));
4027  if (nearest>value) { 
4028    whichWay=1;
4029  } else {
4030    whichWay=0;
4031  }
4032  infeasibility_ = fabs(value-nearest);
4033  bool satisfied=false;
4034  if (infeasibility_<=info->integerTolerance_) {
4035    otherInfeasibility_ = 1.0;
4036    satisfied=true;
4037    if (info->lower_[columnNumber_]!=info->upper_[columnNumber_])
4038      infeasibility_ = 1.0e-5;
4039    else
4040      infeasibility_ = 0.0;
4041  } else if (info->defaultDual_<0.0) {
4042    otherInfeasibility_ = 1.0-infeasibility_;
4043  } else {
4044    const double * pi = info->pi_;
4045    const double * activity = info->rowActivity_;
4046    const double * lower = info->rowLower_;
4047    const double * upper = info->rowUpper_;
4048    const double * element = info->elementByColumn_;
4049    const int * row = info->row_;
4050    const CoinBigIndex * columnStart = info->columnStart_;
4051    const int * columnLength = info->columnLength_;
4052    double direction = info->direction_;
4053    double downMovement = value - floor(value);
4054    double upMovement = 1.0-downMovement;
4055    double valueP = info->objective_[columnNumber_]*direction;
4056    CoinBigIndex start = columnStart[columnNumber_];
4057    CoinBigIndex end = start + columnLength[columnNumber_];
4058    double upEstimate = 0.0;
4059    double downEstimate = 0.0;
4060    if (valueP>0.0)
4061      upEstimate = valueP*upMovement;
4062    else
4063      downEstimate -= valueP*downMovement;
4064    double tolerance = info->primalTolerance_;
4065    for (CoinBigIndex j=start;j<end;j++) {
4066      int iRow = row[j];
4067      if (lower[iRow]<-1.0e20) 
4068        assert (pi[iRow]<=1.0e-3);
4069      if (upper[iRow]>1.0e20) 
4070        assert (pi[iRow]>=-1.0e-3);
4071      valueP = pi[iRow]*direction;
4072      double el2 = element[j];
4073      double value2 = valueP*el2;
4074      double u=0.0;
4075      double d=0.0;
4076      if (value2>0.0)
4077        u = value2;
4078      else
4079        d = -value2;
4080      // if up makes infeasible then make at least default
4081      double newUp = activity[iRow] + upMovement*el2;
4082      if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance)
4083        u = CoinMax(u,info->defaultDual_);
4084      upEstimate += u*upMovement;
4085      // if down makes infeasible then make at least default
4086      double newDown = activity[iRow] - downMovement*el2;
4087      if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance)
4088        d = CoinMax(d,info->defaultDual_);
4089      downEstimate += d*downMovement;
4090    }
4091    if (downEstimate>=upEstimate) {
4092      infeasibility_ = CoinMax(1.0e-12,upEstimate);
4093      otherInfeasibility_ = CoinMax(1.0e-12,downEstimate);
4094      whichWay = 1;
4095    } else {
4096      infeasibility_ = CoinMax(1.0e-12,downEstimate);
4097      otherInfeasibility_ = CoinMax(1.0e-12,upEstimate);
4098      whichWay = 0;
4099    }
4100  }
4101  if (preferredWay_>=0&&!satisfied)
4102    whichWay = preferredWay_;
4103  whichWay_=whichWay;
4104  return infeasibility_;
4105}
4106// Creates a branching object
4107OsiBranchingObject * 
4108OsiSimpleFixedInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const 
4109{
4110  double value = info->solution_[columnNumber_];
4111  value = CoinMax(value, info->lower_[columnNumber_]);
4112  value = CoinMin(value, info->upper_[columnNumber_]);
4113  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
4114  double nearest = floor(value+0.5);
4115  double integerTolerance = info->integerTolerance_;
4116  if (fabs(value-nearest)<integerTolerance) {
4117    // adjust value
4118    if (nearest!=info->upper_[columnNumber_])
4119      value = nearest+2.0*integerTolerance;
4120    else
4121      value = nearest-2.0*integerTolerance;
4122  }
4123  OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way,
4124                                             value);
4125  return branch;
4126}
4127
4128#include <cstdlib>
4129#include <cstdio>
4130#include <cmath>
4131#include <cfloat>
4132#include <cassert>
4133#include <iostream>
4134//#define CGL_DEBUG 2
4135#include "CoinHelperFunctions.hpp"
4136#include "CoinPackedVector.hpp"
4137#include "OsiRowCutDebugger.hpp"
4138#include "CoinWarmStartBasis.hpp"
4139//#include "CglTemporary.hpp"
4140#include "CoinFinite.hpp"
4141//-------------------------------------------------------------------
4142// Generate Stored cuts
4143//-------------------------------------------------------------------
4144void 
4145CglTemporary::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
4146                             const CglTreeInfo info) const
4147{
4148  // Get basic problem information
4149  const double * solution = si.getColSolution();
4150  int numberRowCuts = cuts_.sizeRowCuts();
4151  for (int i=0;i<numberRowCuts;i++) {
4152    const OsiRowCut * rowCutPointer = cuts_.rowCutPtr(i);
4153    double violation = rowCutPointer->violated(solution);
4154    if (violation>=requiredViolation_) 
4155      cs.insert(*rowCutPointer);
4156  }
4157  // delete
4158  cuts_ = OsiCuts();
4159}
4160
4161//-------------------------------------------------------------------
4162// Default Constructor
4163//-------------------------------------------------------------------
4164CglTemporary::CglTemporary ()
4165:
4166CglStored()
4167{
4168}
4169
4170//-------------------------------------------------------------------
4171// Copy constructor
4172//-------------------------------------------------------------------
4173CglTemporary::CglTemporary (const CglTemporary & source) :
4174  CglStored(source)
4175{ 
4176}
4177
4178//-------------------------------------------------------------------
4179// Clone
4180//-------------------------------------------------------------------
4181CglCutGenerator *
4182CglTemporary::clone() const
4183{
4184  return new CglTemporary(*this);
4185}
4186
4187//-------------------------------------------------------------------
4188// Destructor
4189//-------------------------------------------------------------------
4190CglTemporary::~CglTemporary ()
4191{
4192}
4193
4194//----------------------------------------------------------------
4195// Assignment operator
4196//-------------------------------------------------------------------
4197CglTemporary &
4198CglTemporary::operator=(const CglTemporary& rhs)
4199{
4200  if (this != &rhs) {
4201    CglStored::operator=(rhs);
4202  }
4203  return *this;
4204}
4205#endif
Note: See TracBrowser for help on using the repository browser.