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

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

for nonlinear

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