source: trunk/Cbc/src/CbcLinked.cpp

Last change on this file was 2546, checked in by forrest, 4 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
RevLine 
[1854]1/* $Id: CbcLinked.cpp 2546 2019-04-09 13:44:03Z unxusr $ */
[640]2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1573]4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
[640]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
[2464]14static int decodeBit(char *phrase, char *&nextPhrase, double &coefficient, bool ifFirst, const CoinModel &model)
[640]15{
[2464]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;
[640]29    }
[2464]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++;
[1103]37#ifndef NDEBUG
[2464]38      char x = *(pos3 - 1);
39      assert((x >= '0' && x <= '9') || x == '.' || x == '+' || x == '-' || x == 'e');
[1103]40#endif
[640]41    }
42    char saved = *pos2;
[1286]43    *pos2 = '\0';
[2464]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++;
[640]53    }
[2464]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++;
[1103]73#ifndef NDEBUG
[2464]74        char x = *(pos3 - 1);
75        assert((x >= '0' && x <= '9') || x == '.' || x == '+' || x == '-' || x == 'e');
[1103]76#endif
[2464]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();
[640]87    }
[2464]88  }
89  *pos2 = saved;
90  pos = pos2;
91  coefficient = value;
92  nextPhrase = pos;
93  return jColumn;
[640]94}
[687]95#include "ClpQuadraticObjective.hpp"
[640]96#include <cassert>
97#if defined(_MSC_VER)
98// Turn off compiler warning about long names
[2464]99#pragma warning(disable : 4786)
[640]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{
[2464]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_);
[640]142    }
[2464]143    int updated = updateCoefficients(modelPtr_, temp);
[2531]144    //if (updated || 1) {
[2464]145      temp->removeGaps(1.0e-14);
146      ClpMatrixBase *save = modelPtr_->clpMatrix();
[2467]147      ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(save);
[2464]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);
[2531]160    //} else {
161    //  delete temp;
162    //}
[2464]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;
[1286]201        }
[2464]202      }
[640]203    }
[2464]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];
[1286]218        }
[2464]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();
[2546]231        //printf("better qp objective of %g\n", bestObjectiveValue_);
[2464]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();
[2467]240            CglStored *gen2 = dynamic_cast< CglStored * >(gen);
[2464]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;
[1286]257                }
[2464]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;
[1286]265            }
[2464]266          }
267          cbcModel_->unlockThread();
[1286]268        }
[2464]269      }
[640]270    }
[2464]271  }
[640]272}
273//#define WRITE_MATRIX
274#ifdef WRITE_MATRIX
[1286]275static int xxxxxx = 0;
[640]276#endif
277//-----------------------------------------------------------------------------
278void OsiSolverLink::resolve()
279{
[2464]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");
[771]289    }
[2464]290  }
[2546]291  int saveLogLevel=modelPtr_->logLevel();
292  if (saveLogLevel==1)
293    modelPtr_->setLogLevel(0);
[2464]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;
[640]311    }
[2464]312    int updated = updateCoefficients(modelPtr_, temp);
313    if (updated) {
314      temp->removeGaps(1.0e-14);
315      ClpMatrixBase *save = modelPtr_->clpMatrix();
[2467]316      ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(save);
[2464]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  }
[1286]339#ifdef WRITE_MATRIX
[2464]340  {
341    xxxxxx++;
342    char temp[50];
343    sprintf(temp, "bb%d", xxxxxx);
344    writeMps(temp);
345    printf("wrote bb%d\n", xxxxxx);
346  }
[640]347#endif
[2464]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      }
[640]357    }
[2464]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          }
[1286]395        }
[2464]396        returnCode = allFixed ? fathom(allFixed) : 0;
397      } else {
398        returnCode = 0;
399      }
400    } else {
401      returnCode = 0;
[640]402    }
[2464]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;
[1286]420        }
[2464]421      }
[640]422    }
[2464]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
[1103]429#ifndef NDEBUG
[2464]430        double direction = modelPtr_->optimizationDirection();
431        assert(direction == 1.0);
[1103]432#endif
[2464]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++) {
[2467]440          OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]441          if (obj) {
442            value += obj->xyCoefficient(solution);
443          }
444        }
445        if (value < bestObjectiveValue_ - 1.0e-3) {
[2546]446          // printf("obj of %g\n", value);
[2464]447          //modelPtr_->setDualObjectiveLimit(value);
448          delete[] bestSolution_;
449          bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(), modelPtr_->getNumCols());
450          bestObjectiveValue_ = value;
451          if (maxIts <= 10000 && cbcModel_) {
[2467]452            OsiSolverLink *solver2 = dynamic_cast< OsiSolverLink * >(cbcModel_->solver());
[2464]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();
[2467]470              CglStored *gen2 = dynamic_cast< CglStored * >(gen);
[2464]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
[1286]478                int i;
[2464]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++) {
[2467]490                  OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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];
[1286]503                    }
[2464]504                  }
[1286]505                }
[2464]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                  }
[1286]517                }
[2464]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++) {
[2467]535          OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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        }
[2531]547#if 0
548        if (numberContinuous) {
[2464]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++) {
[2467]558              OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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);
[1286]587                    }
[2464]588                    // see if problem
[1103]589#ifndef NDEBUG
[2464]590                    double lambda[4];
[1103]591#endif
[2464]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];
[1103]600#ifndef NDEBUG
[2464]601                    double infeasibility = obj->computeLambdas(xB, yB, xybar, lambda);
602                    assert(infeasibility < 1.0e-9);
[1103]603#endif
[2464]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
[1103]619#ifndef NDEBUG
[2464]620                    double lambda[4];
[1103]621#endif
[2464]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];
[1103]628#ifndef NDEBUG
[2464]629                    double infeasibility = obj->computeLambdas(xB, yB, xybar, lambda);
630                    assert(infeasibility < 1.0e-9);
[1103]631#endif
[2464]632                    setColLower(yColumn, newLower);
633                    setColUpper(yColumn, newUpper);
634                  }
635                  newGap = CoinMax(newGap, CoinMax(gapX, gapY));
[1286]636                }
[2464]637              }
[1286]638            }
[2464]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        }
[2531]650#endif
[2464]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();
[2467]664            CglTemporary *gen2 = dynamic_cast< CglTemporary * >(gen);
[2464]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];
[1286]678                }
[2464]679                qpTemp.setLogLevel(modelPtr_->logLevel());
680                qpTemp.primal();
[2546]681                // assert(!qpTemp.problemStatus());
[2464]682                if (qpTemp.objectiveValue() < bestObjectiveValue_ - 1.0e-3 && !qpTemp.problemStatus()) {
683                  solution2 = CoinCopyOfArray(qpTemp.primalColumnSolution(), numberColumns);
684                } else {
[2546]685                  // printf("QP says expensive - kill\n");
[2464]686                  modelPtr_->setProblemStatus(1);
687                  modelPtr_->setObjectiveValue(COIN_DBL_MAX);
688                  break;
[1286]689                }
[2464]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;
[1286]699            }
[2464]700          }
[1286]701        }
[2464]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();
[2467]709          CglTemporary *gen2 = dynamic_cast< CglTemporary * >(gen);
[2464]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++) {
[2467]738                OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[whichNonLinear_[i]]);
[2464]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      }
[640]787    }
[2464]788  } else {
789    modelPtr_->setProblemStatus(1);
790    modelPtr_->setObjectiveValue(COIN_DBL_MAX);
791  }
[2546]792  modelPtr_->setLogLevel(saveLogLevel);
[640]793}
[700]794// Do OA cuts
[2464]795int OsiSolverLink::doAOCuts(CglTemporary *cutGen, const double *solution, const double *solution2)
[700]796{
[2464]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++) {
[2467]817    OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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      }
[700]831    }
[2464]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;
[700]843    }
[2464]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;
[700]859}
[640]860
861//#############################################################################
862// Constructors, destructors clone and assignment
863//#############################################################################
864
865//-------------------------------------------------------------------
[1286]866// Default Constructor
[640]867//-------------------------------------------------------------------
[2464]868OsiSolverLink::OsiSolverLink()
869  : CbcOsiSolver()
[640]870{
[2464]871  gutsOfDestructor(true);
[640]872}
[1393]873#ifdef JJF_ZERO
[640]874/* returns
875   sequence of nonlinear or
876   -1 numeric
877   -2 not found
878   -3 too many terms
879*/
[2464]880static int getVariable(const CoinModel &model, char *expression,
881  int &linear)
[640]882{
[2464]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;
[1286]897        }
[2464]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);
[1286]915            first = strstr(expression, name);
916            if (first) {
[2464]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              }
[1286]925            }
[2464]926          }
927          break;
[1286]928        }
[2464]929      }
[640]930    }
[2464]931    if (non == -1)
932      non = -2;
933  }
934  return non;
[640]935}
936#endif
[1286]937/* This creates from a coinModel object
[640]938
939   if errors.then number of sets is -1
[1286]940
[640]941   This creates linked ordered sets information.  It assumes -
942
943   for product terms syntax is yy*f(zz)
[1286]944   also just f(zz) is allowed
[640]945   and even a constant
946
947   modelObject not const as may be changed as part of process.
948*/
[2464]949OsiSolverLink::OsiSolverLink(CoinModel &coinModel)
950  : CbcOsiSolver()
[640]951{
[2464]952  gutsOfDestructor(true);
953  load(coinModel);
[640]954}
955// need bounds
[2464]956static void fakeBounds(OsiSolverInterface *solver, int column, double maximumValue,
957  CoinModel *model1, CoinModel *model2)
[640]958{
[2464]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;
[1286]995    }
[2464]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++;
[1286]1004    }
[2464]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;
[1286]1014    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
[2464]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;
[1286]1071        }
[2464]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();
[1286]1104            }
[2464]1105            ifFirst = false;
1106          }
1107          coinModelOriginal.setElement(iRow, iColumn, linearTerm);
[1286]1108        }
[2464]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_);
[1286]1128        }
[2464]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];
[1286]1154            if (coinModel_.isInteger(iColumn))
[2464]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;
[1286]1175        }
[2464]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;
[1286]1209        }
[2464]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;
[2467]1229    const OsiObject **justBi = const_cast< const OsiObject ** >(objects + nInt);
[2464]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
[1286]1254            }
[2464]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]++;
[1286]1274            }
[2464]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;
[1286]1282            }
[2464]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;
[1286]1293        }
[2464]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();
[1286]1361            }
[2464]1362            ifFirst = false;
1363          }
[1286]1364        }
[2464]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++) {
[2467]1427        OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]1428        if (obj) {
1429          int xyRow = obj->xyRow();
1430          assert(xyRow >= 0 && xyRow < numberRows2); // even if obj we should move
1431          whichRows[xyRow]++;
[1286]1432        }
[2464]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;
[1286]1444        }
[2464]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;
[1286]1456        }
[2464]1457      }
1458      whichNonLinear_ = new int[n];
1459      for (i = 0; i < numberObjects_; i++) {
[2467]1460        OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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;
[1286]1467        }
[2464]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++) {
[2467]1484    OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[iObject]);
[2464]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_
[1286]1508        }
[2464]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?
[1286]1523        } else {
[2464]1524          // we would need extra code
1525          abort();
[1286]1526        }
[2464]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      }
[640]1533    }
[2464]1534  }
1535  delete[] which;
1536  if ((specialOptions2_ & 16) != 0)
1537    addTighterConstraints();
[687]1538}
1539// Add reformulated bilinear constraints
[2464]1540void OsiSolverLink::addTighterConstraints()
[687]1541{
[2464]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++) {
[2467]1559    OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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++;
[687]1572    }
[2464]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];
[1103]1596#ifndef NDEBUG
[2464]1597      const double *columnLower = getColLower();
[1103]1598#endif
[2464]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;
[1286]1616        }
[2464]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      }
[687]1674    }
[2464]1675  }
[1393]1676#ifdef JJF_ZERO
[2464]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;
[1286]1692        }
[2464]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      }
[687]1725    }
[2464]1726  }
[640]1727#endif
[2464]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;
[640]1738}
1739// Set all biLinear priorities on x-x variables
[2464]1740void OsiSolverLink::setBiLinearPriorities(int value, double meshSize)
[640]1741{
[2464]1742  OsiObject **newObject = new OsiObject *[numberObjects_];
1743  int numberOdd = 0;
1744  int i;
1745  for (i = 0; i < numberObjects_; i++) {
[2467]1746    OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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      }
[640]1765    }
[2464]1766  }
1767  addObjects(numberOdd, newObject);
1768  for (i = 0; i < numberOdd; i++)
1769    delete newObject[i];
1770  delete[] newObject;
[640]1771}
[687]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*/
[2464]1779void OsiSolverLink::setBranchingStrategyOnVariables(int strategyValue, int priorityValue,
1780  int mode)
[687]1781{
[2464]1782  int i;
1783  for (i = 0; i < numberObjects_; i++) {
[2467]1784    OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]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      }
[687]1801    }
[2464]1802  }
[687]1803}
1804
[640]1805// Say convex (should work it out)
[2464]1806void OsiSolverLink::sayConvex(bool convex)
[1286]1807{
[2464]1808  specialOptions2_ |= 4;
1809  if (convex_) {
1810    for (int iNon = 0; iNon < numberNonLinearRows_; iNon++) {
1811      convex_[iNon] = convex ? 1 : -1;
[640]1812    }
[2464]1813  }
[640]1814}
1815// Set all mesh sizes on x-x variables
[2464]1816void OsiSolverLink::setMeshSizes(double value)
[640]1817{
[2464]1818  int i;
1819  for (i = 0; i < numberObjects_; i++) {
[2467]1820    OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[i]);
[2464]1821    if (obj) {
1822      if (obj->xMeshSize() < 1.0 && obj->yMeshSize() < 1.0) {
[1393]1823#ifdef JJF_ZERO
[2464]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));
[640]1830#endif
[2464]1831        obj->setMeshSizes(this, value, value);
1832      }
[640]1833    }
[2464]1834  }
[640]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*/
[1286]1842double *
1843OsiSolverLink::nonlinearSLP(int numberPasses, double deltaTolerance)
[640]1844{
[2464]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();
[1286]1886        }
[2464]1887        ifFirst = false;
1888      }
[640]1889    }
[2464]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;
[1286]1913        }
[2464]1914      }
1915      triple = coinModel.next(triple);
1916      n++;
[640]1917    }
[2464]1918    if (!linear) {
1919      markNonlinear[iColumn] = 1;
[640]1920    }
[2464]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;
[640]1940    ClpSimplex tempModel;
[1286]1941    tempModel.loadProblem(coinModel, true);
[2464]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++;
[640]1960    }
[2464]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      }
[640]1985    }
[2464]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();
[1286]2011
[2464]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
[1286]2055    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
[2464]2056      iColumn = listNonLinearColumn[jNon];
2057      coinModel.associateElement(coinModel.columnName(iColumn), solution[iColumn]);
[640]2058    }
[2464]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;
[640]2076    }
[2464]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));
[1286]2098            }
[2464]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);
[1286]2105            }
[2464]2106          }
[1286]2107        } else {
[2464]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));
[1286]2115            }
[2464]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);
[1286]2122            }
[2464]2123          }
2124          jNon++;
[1286]2125        }
[2464]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          }
[1286]2149        }
[2464]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);
[1286]2172        }
[2464]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++;
[1286]2240        } else {
[2464]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++;
[1286]2244        }
[2464]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    }
[640]2251#ifdef CLP_DEBUG
[2464]2252    if (logLevel & 32)
2253      std::cout << "largest gap is " << maxGap << " "
2254                << numberSmaller + numberSmaller2 << " reduced ("
2255                << numberSmaller << " badMove ), "
2256                << numberLarger << " increased" << std::endl;
[640]2257#endif
[2464]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    }
[1393]2308#ifdef JJF_ZERO
[2464]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    }
[640]2332#endif
[2464]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);
[1286]2338
[640]2339#ifdef CLP_DEBUG
[2464]2340      if (logLevel & 32)
2341        std::cout << "Pass - " << iPass
2342                  << ", target drop is " << targetDrop
2343                  << std::endl;
[640]2344#endif
[2464]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      }
[640]2351#ifdef CLP_DEBUG
[2464]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      }
[640]2363#endif
[2464]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);
[640]2373#ifdef CLP_DEBUG
[2464]2374      if (model.status()) {
2375        model.writeMps("xx.mps");
2376      }
[640]2377#endif
[2464]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
[640]2405#ifdef CLP_DEBUG
[2464]2406      if (logLevel & 32)
2407        printf("Backtracking\n");
[640]2408#endif
[2464]2409      // load old values
2410      for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
[1286]2411        iColumn = listNonLinearColumn[jNon];
[2464]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++) {
[1286]2426        iColumn = listNonLinearColumn[jNon];
[2464]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;
[1286]2434    }
[2464]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());
[640]2488}
2489/* Solve linearized quadratic objective branch and bound.
2490   Return cutoff and OA cut
2491*/
[1286]2492double
[2464]2493OsiSolverLink::linearizedBAB(CglStored *cut)
[640]2494{
[2464]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();
[2467]2518    OsiClpSolverInterface *osiclpModel = dynamic_cast< OsiClpSolverInterface * >(osiModel);
[2464]2519    ClpSimplex *clpModel = osiclpModel->getModelPtr();
[1286]2520
[2464]2521    // Set changed values
[1286]2522
[640]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);
[1286]2531    cbcModel->addCutGenerator(&probing, -1, "Probing", true, false, false, -100, -1, -1);
[2464]2532    cbcModel->cutGenerator(0)->setTiming(true);
[1286]2533
[640]2534    CglGomory gomory;
2535    gomory.setLimitAtRoot(512);
[1286]2536    cbcModel->addCutGenerator(&gomory, -98, "Gomory", true, false, false, -100, -1, -1);
[2464]2537    cbcModel->cutGenerator(1)->setTiming(true);
[1286]2538
[640]2539    CglKnapsackCover knapsackCover;
[1286]2540    cbcModel->addCutGenerator(&knapsackCover, -98, "KnapsackCover", true, false, false, -100, -1, -1);
[2464]2541    cbcModel->cutGenerator(2)->setTiming(true);
[1286]2542
[640]2543    CglClique clique;
2544    clique.setStarCliqueReport(false);
2545    clique.setRowCliqueReport(false);
2546    clique.setMinViolation(0.1);
[1286]2547    cbcModel->addCutGenerator(&clique, -98, "Clique", true, false, false, -100, -1, -1);
[2464]2548    cbcModel->cutGenerator(3)->setTiming(true);
2549
[640]2550    CglMixedIntegerRounding2 mixedIntegerRounding2;
[1286]2551    cbcModel->addCutGenerator(&mixedIntegerRounding2, -98, "MixedIntegerRounding2", true, false, false, -100, -1, -1);
[2464]2552    cbcModel->cutGenerator(4)->setTiming(true);
[1286]2553
[640]2554    CglFlowCover flowCover;
[1286]2555    cbcModel->addCutGenerator(&flowCover, -98, "FlowCover", true, false, false, -100, -1, -1);
[2464]2556    cbcModel->cutGenerator(5)->setTiming(true);
[1286]2557
[640]2558    CglTwomir twomir;
2559    twomir.setMaxElements(250);
[1286]2560    cbcModel->addCutGenerator(&twomir, -99, "Twomir", true, false, false, -100, -1, -1);
[640]2561    cbcModel->cutGenerator(6)->setTiming(true);
[2464]2562    // For now - switch off most heuristics (because CglPreProcess is bad with QP)
2563#ifndef JJF_ONE
[640]2564    CbcHeuristicFPump heuristicFPump(*cbcModel);
[2464]2565    heuristicFPump.setWhen(13);
[640]2566    heuristicFPump.setMaximumPasses(20);
[2464]2567    heuristicFPump.setMaximumRetries(7);
2568    heuristicFPump.setAbsoluteIncrement(4332.64);
[640]2569    cbcModel->addHeuristic(&heuristicFPump);
[2464]2570    heuristicFPump.setInitialWeight(1);
[1286]2571
[640]2572    CbcHeuristicLocal heuristicLocal(*cbcModel);
2573    heuristicLocal.setSearchType(1);
2574    cbcModel->addHeuristic(&heuristicLocal);
[1286]2575
[640]2576    CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
2577    cbcModel->addHeuristic(&heuristicGreedyCover);
[1286]2578
[640]2579    CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
2580    cbcModel->addHeuristic(&heuristicGreedyEquality);
[2464]2581#endif
[1286]2582
[2464]2583    CbcRounding rounding(*cbcModel);
2584    rounding.setHeuristicName("rounding");
2585    cbcModel->addHeuristic(&rounding);
2586
[640]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();
[2464]2595    clpModel->setDualBound(1.0001e+08);
[640]2596    clpModel->setPerturbation(50);
2597    osiclpModel->setSpecialOptions(193);
2598    osiclpModel->messageHandler()->setLogLevel(0);
[1286]2599    osiclpModel->setIntParam(OsiMaxNumIterationHotStart, 100);
2600    osiclpModel->setHintParam(OsiDoReducePrint, true, OsiHintTry);
[640]2601    // You can save some time by switching off message building
2602    // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
[2464]2603
[640]2604    // Solve
[1286]2605
[640]2606    cbcModel->initialSolve();
[2464]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();
[2467]2615    OsiSolverLinearizedQuadratic *solver3 = dynamic_cast< OsiSolverLinearizedQuadratic * >(model2.solver());
[2464]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++) {
[2467]2673      OsiSimpleInteger *obj = dynamic_cast< OsiSimpleInteger * >(object_[iObject]);
[2464]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++) {
[2467]2700      OsiSimpleInteger *obj = dynamic_cast< OsiSimpleInteger * >(object_[iObject]);
[2464]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      }
[2467]2708      OsiBiLinear *objB = dynamic_cast< OsiBiLinear * >(object_[iObject]);
[2464]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++) {
[2467]2731      OsiSimpleInteger *obj = dynamic_cast< OsiSimpleInteger * >(object_[iObject]);
[2464]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      }
[2467]2739      OsiBiLinear *objB = dynamic_cast< OsiBiLinear * >(object_[iObject]);
[2464]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();
[2467]2760  OsiClpSolverInterface *osiclpModel = dynamic_cast< OsiClpSolverInterface * >(osiModel);
[2464]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);
[1886]2844#ifndef NDEBUG
[2464]2845  int returnCode = 0;
[1886]2846#endif
[2464]2847  if (clpModel->tightenPrimalBounds() != 0) {
2848    clpModel->setLogLevel(saveLogLevel);
[1886]2849#ifndef NDEBUG
[2464]2850    returnCode = -1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
[1886]2851#endif
[2464]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;
[1886]2890#ifndef NDEBUG
[2464]2891      returnCode = -1;
[1886]2892#endif
[2464]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      }
[640]2917    }
[2464]2918  }
2919  assert(!returnCode);
2920  abort();
2921  return solution;
[640]2922}
2923// Analyze constraints to see which are convex (quadratic)
[2464]2924void OsiSolverLink::analyzeObjects()
[640]2925{
[2464]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++) {
[2467]2941      OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[whichNonLinear_[i]]);
[2464]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++;
[1286]2971        }
[2464]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          }
[1286]2994        }
[2464]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;
[1286]3015                }
[2464]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              }
[1286]3028            }
[2464]3029          } else {
3030            numberLong++;
3031          }
[1286]3032        }
[2464]3033      }
[640]3034    }
[2464]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;
[640]3063}
3064//-------------------------------------------------------------------
3065// Clone
3066//-------------------------------------------------------------------
[1286]3067OsiSolverInterface *
[1271]3068OsiSolverLink::clone(bool /*copyData*/) const
[640]3069{
[2464]3070  //assert (copyData);
3071  OsiSolverLink *newModel = new OsiSolverLink(*this);
3072  return newModel;
[640]3073}
3074
3075//-------------------------------------------------------------------
[1286]3076// Copy constructor
[640]3077//-------------------------------------------------------------------
[2464]3078OsiSolverLink::OsiSolverLink(
3079  const OsiSolverLink &rhs)
3080  : OsiSolverInterface(rhs)
3081  , CbcOsiSolver(rhs)
[640]3082{
[2464]3083  gutsOfDestructor(true);
3084  gutsOfCopy(rhs);
3085  // something odd happens - try this
3086  OsiSolverInterface::operator=(rhs);
[640]3087}
3088
3089//-------------------------------------------------------------------
[1286]3090// Destructor
[640]3091//-------------------------------------------------------------------
[2464]3092OsiSolverLink::~OsiSolverLink()
[640]3093{
[2464]3094  gutsOfDestructor();
[640]3095}
3096
3097//-------------------------------------------------------------------
[1286]3098// Assignment operator
[640]3099//-------------------------------------------------------------------
3100OsiSolverLink &
[2464]3101OsiSolverLink::operator=(const OsiSolverLink &rhs)
[640]3102{
[2464]3103  if (this != &rhs) {
3104    gutsOfDestructor();
3105    CbcOsiSolver::operator=(rhs);
3106    gutsOfCopy(rhs);
3107  }
3108  return *this;
[640]3109}
[2464]3110void OsiSolverLink::gutsOfDestructor(bool justNullify)
[640]3111{
[2464]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;
[640]3145}
[2464]3146void OsiSolverLink::gutsOfCopy(const OsiSolverLink &rhs)
[640]3147{
[2464]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]);
[640]3172    }
[2464]3173    if (rhs.bestSolution_) {
3174      bestSolution_ = CoinCopyOfArray(rhs.bestSolution_, modelPtr_->getNumCols());
[640]3175    } else {
[2464]3176      bestSolution_ = NULL;
[640]3177    }
[2464]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_);
[640]3192}
3193// Add a bound modifier
[2464]3194void OsiSolverLink::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, int whichVariableAffected,
3195  double multiplier)
[640]3196{
[2464]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;
[640]3203    }
[2464]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);
[640]3215}
3216// Update coefficients
[2464]3217int OsiSolverLink::updateCoefficients(ClpSimplex *solver, CoinPackedMatrix *matrix)
[640]3218{
[2464]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++) {
[2467]3224    OsiBiLinear *obj = dynamic_cast< OsiBiLinear * >(object_[iObject]);
[2464]3225    if (obj) {
3226      numberChanged += obj->updateCoefficients(lower, upper, objective, matrix, &basis_);
[640]3227    }
[2464]3228  }
3229  return numberChanged;
[640]3230}
3231// Set best solution found internally
[2464]3232void OsiSolverLink::setBestSolution(const double *solution, int numberColumns)
[640]3233{
[2464]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));
[640]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*/
[2464]3243void OsiSolverLink::setFixedPriority(int priorityValue)
[640]3244{
[2464]3245  delete[] fixVariables_;
3246  fixVariables_ = NULL;
3247  numberFix_ = 0;
3248  int i;
3249  for (i = 0; i < numberObjects_; i++) {
[2467]3250    OsiSimpleInteger *obj = dynamic_cast< OsiSimpleInteger * >(object_[i]);
[2464]3251    if (obj) {
[1103]3252#ifndef NDEBUG
[2464]3253      int iColumn = obj->columnNumber();
3254      assert(iColumn >= 0);
[1103]3255#endif
[2464]3256      if (obj->priority() < priorityValue)
3257        numberFix_++;
[640]3258    }
[2464]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++) {
[2467]3269      OsiSimpleInteger *obj = dynamic_cast< OsiSimpleInteger * >(object_[i]);
[2464]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          }
[1286]3280        }
[2464]3281      }
[640]3282    }
[2464]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  }
[640]3295}
3296// Gets correct form for a quadratic row - user to delete
[1286]3297CoinPackedMatrix *
[2464]3298OsiSolverLink::quadraticRow(int rowNumber, double *linearRow) const
[640]3299{
[2464]3300  int numberColumns = coinModel_.numberColumns();
3301  CoinZeroN(linearRow, numberColumns);
3302  int numberElements = 0;
[1103]3303#ifndef NDEBUG
[2464]3304  int numberRows = coinModel_.numberRows();
3305  assert(rowNumber >= 0 && rowNumber < numberRows);
[1103]3306#endif
[2464]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++;