source: trunk/Cbc/src/CbcLinked.cpp

Last change on this file was 2546, checked in by forrest, 3 months ago

try and improve qp

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