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

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

nonlinear stuff

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