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

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

quadratic

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