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

Last change on this file since 520 was 520, checked in by forrest, 12 years ago

fix bug in ampl interface and try slp

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