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

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

for nonlinear

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