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

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

for tighter linear constraints in bilinear

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