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

Last change on this file since 1424 was 1393, checked in by lou, 10 years ago

Mark #if 0 with JJF_ZERO and #if 1 with JJF_ONE as a historical reference
point.

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