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

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

small changes for LOS

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