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

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

obviously needs more work

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