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

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

for Jon

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