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

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

add time and try two level bilinear

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