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

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

probably breaking everything

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