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

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

mostly for rins

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