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

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

for Matteo and a bit for Jon

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