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

Last change on this file since 2081 was 1899, checked in by stefan, 7 years ago

fixup svn properties

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 311.5 KB
RevLine 
[1854]1/* $Id: CbcLinked.cpp 1899 2013-04-09 18:12:08Z forrest $ */
[640]2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1573]4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
[640]6#include "CbcConfig.h"
7
8#include "CoinTime.hpp"
9
10#include "CoinHelperFunctions.hpp"
11#include "CoinModel.hpp"
12#include "ClpSimplex.hpp"
13// returns jColumn (-2 if linear term, -1 if unknown) and coefficient
14static
15int decodeBit(char * phrase, char * & nextPhrase, double & coefficient, bool ifFirst, const CoinModel & model)
16{
[1286]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++;
[640]32    }
[1286]33    // if * must be number otherwise must be name
34    if (*pos2 == '*') {
35        char * pos3 = pos;
36        while (pos3 != pos2) {
37            pos3++;
[1103]38#ifndef NDEBUG
[1286]39            char x = *(pos3 - 1);
40            assert ((x >= '0' && x <= '9') || x == '.' || x == '+' || x == '-' || x == 'e');
[1103]41#endif
[1286]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        }
[640]55    }
56    char saved = *pos2;
[1286]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;
[640]66    }
[1286]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++;
[1103]74#ifndef NDEBUG
[1286]75                char x = *(pos3 - 1);
76                assert ((x >= '0' && x <= '9') || x == '.' || x == '+' || x == '-' || x == 'e');
[1103]77#endif
[1286]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        }
[640]89    }
[1286]90    *pos2 = saved;
91    pos = pos2;
92    coefficient = value;
93    nextPhrase = pos;
94    return jColumn;
[640]95}
[687]96#include "ClpQuadraticObjective.hpp"
[640]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{
[1286]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        }
[640]164    }
[1286]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");
[640]176
177    }
[1286]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        }
[640]273    }
274}
275//#define WRITE_MATRIX
276#ifdef WRITE_MATRIX
[1286]277static int xxxxxx = 0;
[640]278#endif
279//-----------------------------------------------------------------------------
280void OsiSolverLink::resolve()
281{
[1286]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        }
[771]292    }
[1286]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);
[1876]332            modelPtr_->setNewRowCopy(NULL);
333            modelPtr_->setClpScaledMatrix(NULL);
[1286]334        } else {
335            delete temp;
336        }
[640]337    }
[1286]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);
[640]345    }
346#endif
[1286]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");
[640]359    }
[1286]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        }
[640]371    }
[1286]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        }
[640]402    }
[1286]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
[1103]428#ifndef NDEBUG
[1286]429                double direction = modelPtr_->optimizationDirection();
430                assert (direction == 1.0);
[1103]431#endif
[1286]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
[1103]587#ifndef NDEBUG
[1286]588                                        double lambda[4];
[1103]589#endif
[1286]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];
[1103]598#ifndef NDEBUG
[1286]599                                        double infeasibility = obj->computeLambdas(xB, yB, xybar, lambda);
600                                        assert (infeasibility < 1.0e-9);
[1103]601#endif
[1286]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
[1103]617#ifndef NDEBUG
[1286]618                                        double lambda[4];
[1103]619#endif
[1286]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];
[1103]626#ifndef NDEBUG
[1286]627                                        double infeasibility = obj->computeLambdas(xB, yB, xybar, lambda);
628                                        assert (infeasibility < 1.0e-9);
[1103]629#endif
[1286]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);
[640]788    }
789}
[700]790// Do OA cuts
[1286]791int
[700]792OsiSolverLink::doAOCuts(CglTemporary * cutGen, const double * solution, const double * solution2)
793{
[1286]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        }
[700]829    }
[1286]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        }
[700]841    }
[1286]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;
[700]856}
[640]857
858//#############################################################################
859// Constructors, destructors clone and assignment
860//#############################################################################
861
862//-------------------------------------------------------------------
[1286]863// Default Constructor
[640]864//-------------------------------------------------------------------
865OsiSolverLink::OsiSolverLink ()
[1286]866        : CbcOsiSolver()
[640]867{
[1286]868    gutsOfDestructor(true);
[640]869}
[1393]870#ifdef JJF_ZERO
[640]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,
[1286]878                       int & linear)
[640]879{
[1286]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;
[640]930    }
[1286]931    return non;
[640]932}
933#endif
[1286]934/* This creates from a coinModel object
[640]935
936   if errors.then number of sets is -1
[1286]937
[640]938   This creates linked ordered sets information.  It assumes -
939
940   for product terms syntax is yy*f(zz)
[1286]941   also just f(zz) is allowed
[640]942   and even a constant
943
944   modelObject not const as may be changed as part of process.
945*/
946OsiSolverLink::OsiSolverLink ( CoinModel & coinModel)
[1286]947        : CbcOsiSolver()
[640]948{
[1286]949    gutsOfDestructor(true);
950    load(coinModel);
[640]951}
952// need bounds
[1286]953static void fakeBounds(OsiSolverInterface * solver, int column, double maximumValue,
954                       CoinModel * model1, CoinModel * model2)
[640]955{
[1286]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    }
[640]972}
[1286]973void OsiSolverLink::load ( CoinModel & coinModelOriginal, bool tightenBounds, int logLevel)
[640]974{
[1286]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        }
[640]1005    }
[1286]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        }
[640]1470    }
[1286]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        }
[640]1531    }
1532    delete [] which;
[1286]1533    if ((specialOptions2_&16) != 0)
1534        addTighterConstraints();
[687]1535}
1536// Add reformulated bilinear constraints
[1286]1537void
[687]1538OsiSolverLink::addTighterConstraints()
1539{
[1286]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        }
[687]1571    }
[1286]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];
[1103]1594#ifndef NDEBUG
[1286]1595            const double * columnLower = getColLower();
[1103]1596#endif
[1286]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        }
[687]1674    }
[1393]1675#ifdef JJF_ZERO
[1286]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        }
[687]1726    }
[640]1727#endif
[1286]1728    delete [] xW;
1729    delete [] yW;
1730    delete [] wW;
1731    delete [] alphaW;
1732    delete [] addColumn;
1733    delete [] addElement;
1734    delete [] addW;
1735    delete [] mark;
1736    delete [] list;
1737    delete [] objW;
[640]1738}
1739// Set all biLinear priorities on x-x variables
[1286]1740void
1741OsiSolverLink::setBiLinearPriorities(int value, double meshSize)
[640]1742{
[1286]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        }
[640]1767    }
[1286]1768    addObjects(numberOdd, newObject);
1769    for (i = 0; i < numberOdd; i++)
1770        delete newObject[i];
1771    delete [] newObject;
[640]1772}
[687]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*/
[1286]1780void
[687]1781OsiSolverLink::setBranchingStrategyOnVariables(int strategyValue, int priorityValue,
[1286]1782        int mode)
[687]1783{
[1286]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        }
[687]1805    }
1806}
1807
[640]1808// Say convex (should work it out)
[1286]1809void
[640]1810OsiSolverLink::sayConvex(bool convex)
[1286]1811{
1812    specialOptions2_ |= 4;
1813    if (convex_) {
1814        for (int iNon = 0; iNon < numberNonLinearRows_; iNon++) {
1815            convex_[iNon] = convex ? 1 : -1;
1816        }
[640]1817    }
1818}
1819// Set all mesh sizes on x-x variables
[1286]1820void
[640]1821OsiSolverLink::setMeshSizes(double value)
1822{
[1286]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) {
[1393]1828#ifdef JJF_ZERO
[1286]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));
[640]1835#endif
[1286]1836                obj->setMeshSizes(this, value, value);
1837            }
1838        }
[640]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*/
[1286]1847double *
1848OsiSolverLink::nonlinearSLP(int numberPasses, double deltaTolerance)
[640]1849{
[1286]1850    if (!coinModel_.numberRows()) {
1851        printf("Model not set up or nonlinear arrays not created!\n");
1852        return NULL;
[640]1853    }
[1286]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;
[1886]1888                } else if (jColumn != -2) {
[1286]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;
[1886]1913                    } else if (jColumn != -2) {
[1286]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        }
[640]1926    }
[1286]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        }
[640]1939    }
[1286]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;
[640]1950    }
[1286]1951    // Create artificials
[640]1952    ClpSimplex tempModel;
[1286]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        }
[640]1966    }
[1286]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        }
[640]1991    }
[1286]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];
[640]2009    }
[1286]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];
[640]2025    }
[1286]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;
[640]2042    }
2043    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
[1286]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        }
[640]2257#ifdef CLP_DEBUG
[1286]2258        if (logLevel&32)
2259            std::cout << "largest gap is " << maxGap << " "
2260                      << numberSmaller + numberSmaller2 << " reduced ("
2261                      << numberSmaller << " badMove ), "
2262                      << numberLarger << " increased" << std::endl;
[640]2263#endif
[1286]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        }
[1393]2314#ifdef JJF_ZERO
[1286]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        }
[640]2338#endif
[1286]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
[640]2345#ifdef CLP_DEBUG
[1286]2346            if (logLevel&32)
2347                std::cout << "Pass - " << iPass
2348                          << ", target drop is " << targetDrop
2349                          << std::endl;
[640]2350#endif
[1286]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            }
[640]2357#ifdef CLP_DEBUG
[1286]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            }
[640]2369#endif
[1286]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);
[640]2379#ifdef CLP_DEBUG
[1286]2380            if (model.status()) {
2381                model.writeMps("xx.mps");
2382            }
[640]2383#endif
[1286]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
[640]2411#ifdef CLP_DEBUG
[1286]2412            if (logLevel&32)
2413                printf("Backtracking\n");
[640]2414#endif
[1286]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        }
[640]2441    }
[1286]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());
[640]2494}
2495/* Solve linearized quadratic objective branch and bound.
2496   Return cutoff and OA cut
2497*/
[1286]2498double
2499OsiSolverLink::linearizedBAB(CglStored * cut)
[640]2500{
[1286]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)
[1393]2569#ifndef JJF_ONE
[1286]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);
[640]2657    }
[1286]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);
[640]2763    // Now do requested saves and modifications
[1286]2764    CbcModel * cbcModel = & model;
2765    OsiSolverInterface * osiModel = model.solver();
[640]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);
[1286]2775    probing.setRowCuts(0);
[640]2776    probing.setUsingObjective(true);
[1286]2777    cbcModel->addCutGenerator(&probing, -1, "Probing", true, false, false, -100, -1, -1);
2778
[640]2779    CglGomory gomory;
2780    gomory.setLimitAtRoot(512);
[1286]2781    cbcModel->addCutGenerator(&gomory, -98, "Gomory", true, false, false, -100, -1, -1);
2782
[640]2783    CglKnapsackCover knapsackCover;
[1286]2784    cbcModel->addCutGenerator(&knapsackCover, -98, "KnapsackCover", true, false, false, -100, -1, -1);
2785
[640]2786    CglClique clique;
2787    clique.setStarCliqueReport(false);
2788    clique.setRowCliqueReport(false);
2789    clique.setMinViolation(0.1);
[1286]2790    cbcModel->addCutGenerator(&clique, -98, "Clique", true, false, false, -100, -1, -1);
[640]2791    CglMixedIntegerRounding2 mixedIntegerRounding2;
[1286]2792    cbcModel->addCutGenerator(&mixedIntegerRounding2, -98, "MixedIntegerRounding2", true, false, false, -100, -1, -1);
2793
[640]2794    CglFlowCover flowCover;
[1286]2795    cbcModel->addCutGenerator(&flowCover, -98, "FlowCover", true, false, false, -100, -1, -1);
2796
[640]2797    CglTwomir twomir;
2798    twomir.setMaxElements(250);
[1286]2799    cbcModel->addCutGenerator(&twomir, -99, "Twomir", true, false, false, -100, -1, -1);
[640]2800    cbcModel->cutGenerator(6)->setTiming(true);
[1286]2801
[640]2802    CbcHeuristicFPump heuristicFPump(*cbcModel);
[1286]2803    heuristicFPump.setWhen(1);
[640]2804    heuristicFPump.setMaximumPasses(20);
[1286]2805    heuristicFPump.setDefaultRounding(0.5);
[640]2806    cbcModel->addHeuristic(&heuristicFPump);
[1286]2807
2808    CbcRounding rounding(*cbcModel);
2809    cbcModel->addHeuristic(&rounding);
2810
[640]2811    CbcHeuristicLocal heuristicLocal(*cbcModel);
2812    heuristicLocal.setSearchType(1);
2813    cbcModel->addHeuristic(&heuristicLocal);
[1286]2814
[640]2815    CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
2816    cbcModel->addHeuristic(&heuristicGreedyCover);
[1286]2817
[640]2818    CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
2819    cbcModel->addHeuristic(&heuristicGreedyEquality);
[1286]2820
2821    CbcCompareDefault compare;
2822    cbcModel->setNodeComparison(compare);
[640]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);
[1286]2829    clpModel->setNumberIterations(1);
[640]2830    // For branchAndBound this may help
2831    clpModel->defaultFactorizationFrequency();
[1286]2832    clpModel->setDualBound(6.71523e+07);
[640]2833    clpModel->setPerturbation(50);
2834    osiclpModel->setSpecialOptions(193);
2835    osiclpModel->messageHandler()->setLogLevel(0);
[1286]2836    osiclpModel->setIntParam(OsiMaxNumIterationHotStart, 100);
2837    osiclpModel->setHintParam(OsiDoReducePrint, true, OsiHintTry);
[640]2838    // You can save some time by switching off message building
2839    // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
2840    // Solve
[1286]2841
[640]2842    cbcModel->initialSolve();
[1286]2843    //double cutoff = model_->getCutoff();
2844    if (!cbcModel_)
2845        cbcModel->setCutoff(1.0e50);
[640]2846    else
[1286]2847        cbcModel->setCutoff(cbcModel_->getCutoff());
2848    int saveLogLevel = clpModel->logLevel();
2849    clpModel->setLogLevel(0);
[1886]2850#ifndef NDEBUG
[1286]2851    int returnCode = 0;
[1886]2852#endif
[1286]2853    if (clpModel->tightenPrimalBounds() != 0) {
2854        clpModel->setLogLevel(saveLogLevel);
[1886]2855#ifndef NDEBUG
[1286]2856        returnCode = -1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
[1886]2857#endif
[1286]2858        //clpModel->writeMps("infeas2.mps");
[640]2859    } else {
[1286]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;
[1886]2897#ifndef NDEBUG
[1286]2898            returnCode = -1;
[1886]2899#endif
[1286]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        }
[640]2925    }
[1286]2926    assert (!returnCode);
2927    abort();
2928    return solution;
[640]2929}
2930// Analyze constraints to see which are convex (quadratic)
[1286]2931void
[640]2932OsiSolverLink::analyzeObjects()
2933{
[1286]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;
[640]3068    }
[1286]3069    delete [] start;
[640]3070}
3071//-------------------------------------------------------------------
3072// Clone
3073//-------------------------------------------------------------------
[1286]3074OsiSolverInterface *
[1271]3075OsiSolverLink::clone(bool /*copyData*/) const
[640]3076{
[1286]3077    //assert (copyData);
3078    OsiSolverLink * newModel = new OsiSolverLink(*this);
3079    return newModel;
[640]3080}
3081
3082
3083//-------------------------------------------------------------------
[1286]3084// Copy constructor
[640]3085//-------------------------------------------------------------------
3086OsiSolverLink::OsiSolverLink (
[1286]3087    const OsiSolverLink & rhs)
3088        : OsiSolverInterface(rhs),
3089        CbcOsiSolver(rhs)
[640]3090{
[1286]3091    gutsOfDestructor(true);
3092    gutsOfCopy(rhs);
3093    // something odd happens - try this
3094    OsiSolverInterface::operator=(rhs);
[640]3095}
3096
3097//-------------------------------------------------------------------
[1286]3098// Destructor
[640]3099//-------------------------------------------------------------------
3100OsiSolverLink::~OsiSolverLink ()
3101{
[1286]3102    gutsOfDestructor();
[640]3103}
3104
3105//-------------------------------------------------------------------
[1286]3106// Assignment operator
[640]3107//-------------------------------------------------------------------
3108OsiSolverLink &
[1286]3109OsiSolverLink::operator=(const OsiSolverLink & rhs)
[640]3110{
[1286]3111    if (this != &rhs) {
3112        gutsOfDestructor();
3113        CbcOsiSolver::operator=(rhs);
3114        gutsOfCopy(rhs);
3115    }
3116    return *this;
[640]3117}
[1286]3118void
[640]3119OsiSolverLink::gutsOfDestructor(bool justNullify)
3120{
[1286]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;
[640]3154}
[1286]3155void
[640]3156OsiSolverLink::gutsOfCopy(const OsiSolverLink & rhs)
3157{
[1286]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        }
[640]3188    }
[1286]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_);
[640]3198    } else {
[1286]3199        quadraticModel_ = NULL;
[640]3200    }
[1286]3201    fixVariables_ = CoinCopyOfArray(rhs.fixVariables_, numberFix_);
[640]3202}
3203// Add a bound modifier
[1286]3204void
3205OsiSolverLink::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, int whichVariableAffected,
3206                                double multiplier)
[640]3207{
[1286]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        }
[640]3215    }
[1286]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);
[640]3226}
3227// Update coefficients
3228int
3229OsiSolverLink::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
3230{
[1286]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        }
[640]3240    }
[1286]3241    return numberChanged;
[640]3242}
3243// Set best solution found internally
[1286]3244void
[640]3245OsiSolverLink::setBestSolution(const double * solution, int numberColumns)
3246{
[1286]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));
[640]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*/
[1286]3256void
[640]3257OsiSolverLink::setFixedPriority(int priorityValue)
3258{
[1286]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) {
[1103]3266#ifndef NDEBUG
[1286]3267            int iColumn = obj->columnNumber();
3268            assert (iColumn >= 0);
[1103]3269#endif
[1286]3270            if (obj->priority