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

Last change on this file since 1886 was 1886, checked in by stefan, 6 years ago

fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables

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