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

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

update branches/devel for threads

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