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

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

heuristics and nonlinear

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