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

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

nonlinear

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