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

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

changes

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