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

Last change on this file since 1733 was 1573, checked in by lou, 9 years ago

Change to EPL license notice.

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