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

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

for dubiois optionas

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