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

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

trying to move ampl stuff away

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