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

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

changes for gcc 4.3

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