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

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

changes so OsiIntegers? have pseudocosts updated on branch

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