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

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

try and make a bit faster

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