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

Last change on this file since 1876 was 1876, checked in by forrest, 7 years ago

changes for cuts and extra variables

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