source: trunk/Cbc/src/CbcLinked.cpp @ 700

Last change on this file since 700 was 700, checked in by forrest, 14 years ago

changes for postprocess loop and LOS

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