source: trunk/Alps/examples/Abc/AbcModel.cpp @ 277

Last change on this file since 277 was 277, checked in by andreasw, 13 years ago

first working version with autotools

File size: 34.0 KB
Line 
1/*===========================================================================*
2 * This file is part of the Abstract Library for Parallel Search (ALPS).     *
3 *                                                                           *
4 * ALPS is distributed under the Common Public License as part of the        *
5 * COIN-OR repository (http://www.coin-or.org).                              *
6 *                                                                           *
7 * Authors: Yan Xu, SAS Institute Inc.                                       *
8 *          Ted Ralphs, Lehigh University                                    *
9 *          Laszlo Ladanyi, IBM T.J. Watson Research Center                  *
10 *          Matthew Saltzman, Clemson University                             *
11 *                                                                           *
12 *                                                                           *
13 * Copyright (C) 2001-2004, International Business Machines                  *
14 * Corporation, Lehigh University, Yan Xu, Ted Ralphs, Matthew Salzman and   *
15 * others. All Rights Reserved.                                              *
16 *===========================================================================*/
17
18//#############################################################################
19// This file is modified from SbbModel.cpp
20//#############################################################################
21
22#include <iostream>
23
24#include "OsiRowCut.hpp"
25#include "OsiColCut.hpp"
26#include "OsiRowCutDebugger.hpp"
27
28#include "AlpsTreeNode.h"
29
30#include "AbcCutGenerator.h"
31#include "AbcHeuristic.h"
32#include "AbcMessage.h"
33#include "AbcModel.h"
34#include "AbcTreeNode.h"
35#include "AbcNodeDesc.h"
36
37
38//#############################################################################
39//#############################################################################
40
41void
42AbcModel::initialSolve()
43{
44    assert (solver_);
45    solver_->initialSolve();
46}
47
48//#############################################################################
49//  Parameters:
50//  cuts:        (o) all cuts generated in this round of cut generation
51//  numberTries: (i) the maximum number of iterations for this round of cut
52//                   generation; no a priori limit if 0
53//  whichGenerator: (i/o) whichGenerator[i] is loaded with the index of the
54//                        generator that produced cuts[i]; reallocated as
55//                        required
56//  numberOldActiveCuts: (o) the number of active cuts at this node from
57//                           previous rounds of cut generation
58//  numberNewCuts:       (o) the number of cuts produced in this round of cut
59//                           generation
60//  maximumWhich:      (i/o) capacity of whichGenerator; may be updated if
61//                           whichGenerator grows.
62//  cutDuringRampup:    (i) Whether generating cuts during rampup
63//  found: (o)  great than 0 means that heuristics found solutions;
64//              otherwise not.
65bool 
66AbcModel::solveWithCuts(OsiCuts & cuts, int numberTries, 
67                        AbcTreeNode * node, int & numberOldActiveCuts, 
68                        int & numberNewCuts, int & maximumWhich, 
69                        int *& whichGenerator, const bool cutDuringRampup,
70                        int & found)
71{
72    found = -10;
73    bool feasible;
74    int lastNumberCuts = 0;
75    double lastObjective = -1.0e100 ;
76    int violated = 0;
77    int numberRowsAtStart = solver_->getNumRows();
78    int numberColumns = solver_->getNumCols();
79
80    numberOldActiveCuts = numberRowsAtStart - numberRowsAtContinuous_;
81    numberNewCuts = 0;
82   
83    feasible = resolve(); 
84    if(!feasible) {
85        return false;  // If lost feasibility, bail out right now
86    }
87   
88    reducedCostFix();
89    const double *lower = solver_->getColLower();
90    const double *upper = solver_->getColUpper();
91    const double *solution = solver_->getColSolution();
92
93    double minimumDrop = minimumDrop_;
94    if (numberTries < 0) { 
95        numberTries = -numberTries;
96        minimumDrop = -1.0; 
97    }
98
99    //-------------------------------------------------------------------------
100    // Is it time to scan the cuts in order to remove redundant cuts? If so,
101    // set up to do it.
102# define SCANCUTS 100 
103    int *countColumnCuts = NULL;
104    int *countRowCuts = NULL;
105    bool fullScan = false;
106    if ((numberNodes_ % SCANCUTS) == 0) {
107        fullScan = true;
108        countColumnCuts = new int[numberCutGenerators_ + numberHeuristics_];
109        countRowCuts = new int[numberCutGenerators_ + numberHeuristics_];
110        memset(countColumnCuts, 0,
111               (numberCutGenerators_ + numberHeuristics_) * sizeof(int));
112        memset(countRowCuts, 0,
113               (numberCutGenerators_ + numberHeuristics_) * sizeof(int));
114    }
115
116    double direction = solver_->getObjSense();
117    double startObjective = solver_->getObjValue() * direction;
118
119    int numberPasses = 0;
120    double primalTolerance = 1.0e-7;
121
122    //-------------------------------------------------------------------------
123    // Start cut generation loop
124    do {
125        numberPasses++;
126        numberTries--;
127        OsiCuts theseCuts;
128   
129        // First check if there are cuts violated in global cut pool
130        if (numberPasses == 1 && howOftenGlobalScan_ > 0 &&
131            (numberNodes_ % howOftenGlobalScan_) == 0) { 
132            int numberCuts = globalCuts_.sizeColCuts();
133            int i;
134            for ( i = 0; i < numberCuts; ++i) { 
135                const OsiColCut *thisCut = globalCuts_.colCutPtr(i);
136                if (thisCut->violated(solution) > primalTolerance) {
137                    printf("Global cut added - violation %g\n",
138                           thisCut->violated(solution));
139                    theseCuts.insert(*thisCut);
140                }
141            }
142            numberCuts = globalCuts_.sizeRowCuts();
143            for ( i = 0; i < numberCuts; ++i) {
144                const OsiRowCut * thisCut = globalCuts_.rowCutPtr(i);
145                if (thisCut->violated(solution) > primalTolerance) {
146                    printf("Global cut added - violation %g\n",
147                           thisCut->violated(solution));
148                    theseCuts.insert(*thisCut);
149                }
150            }
151        }
152
153        //---------------------------------------------------------------------
154        // Generate new cuts (global and/or local) and/or apply heuristics
155        // NOTE: Make sure CglProbing is added FIRST
156        double * newSolution = new double [numberColumns];
157        double heuristicValue = getCutoff();
158
159#if defined(ABC_DEBUG_MORE)
160            std::cout << "numberCutGenerators_ = " << numberCutGenerators_
161                      << "numberHeuristics_ = " << numberHeuristics_
162                      << std::endl;
163#endif
164        for (int i = 0; i < numberCutGenerators_ + numberHeuristics_; ++i) {
165            int numberRowCutsBefore = theseCuts.sizeRowCuts();
166            int numberColumnCutsBefore = theseCuts.sizeColCuts();
167            if (i < numberCutGenerators_) {
168                if (cutDuringRampup) {
169                    bool mustResolve = 
170                        generator_[i]->generateCuts(theseCuts, fullScan);
171                    if (mustResolve) {
172                        feasible = resolve();
173                        if (!feasible)
174                            break;
175                    }
176                }
177            } 
178            else { 
179                double saveValue = heuristicValue;
180                int ifSol = heuristic_[i-numberCutGenerators_]->
181                    solution(heuristicValue, newSolution);
182                    //    solution(heuristicValue, newSolution, theseCuts);
183
184                if (ifSol > 0) {
185                    found = i;
186                } 
187                else if (ifSol < 0) {
188                    heuristicValue = saveValue;
189                }
190            }
191            int numberRowCutsAfter = theseCuts.sizeRowCuts();
192            int numberColumnCutsAfter = theseCuts.sizeColCuts();
193            int numberBefore =
194                numberRowCutsBefore + numberColumnCutsBefore + lastNumberCuts;
195            int numberAfter =
196                numberRowCutsAfter + numberColumnCutsAfter + lastNumberCuts;
197            if (numberAfter > maximumWhich) {
198                maximumWhich = max(maximumWhich * 2 + 100, numberAfter);
199                int * temp = new int[2 * maximumWhich];
200                memcpy(temp, whichGenerator, numberBefore * sizeof(int));
201                delete [] whichGenerator;
202                whichGenerator = temp;
203            }
204            int j;
205            if (fullScan) {
206                countRowCuts[i] += numberRowCutsAfter - 
207                    numberRowCutsBefore;
208                countColumnCuts[i] += numberColumnCutsAfter - 
209                    numberColumnCutsBefore;
210            }
211            for (j = numberRowCutsBefore; j < numberRowCutsAfter; ++j) {
212                whichGenerator[numberBefore++] = i;
213                const OsiRowCut * thisCut = theseCuts.rowCutPtr(j);
214                if (thisCut->globallyValid()) {
215                    globalCuts_.insert(*thisCut);
216                }
217            }
218            for (j = numberColumnCutsBefore; j < numberColumnCutsAfter; ++j) {
219                whichGenerator[numberBefore++] = i;
220                const OsiColCut * thisCut = theseCuts.colCutPtr(j);
221                if (thisCut->globallyValid()) {
222                    globalCuts_.insert(*thisCut);
223                }
224            }
225        }
226
227        //---------------------------------------------------------------------
228        // If found a solution, Record it before we free the vector
229        if (found >= 0) {
230            bool better = 
231                setBestSolution(ABC_ROUNDING, heuristicValue, newSolution);
232            //    if (!better){
233            //  found = -1;
234            //}
235            //std::cout << "better = "  << better
236            //        << "; found = " << found << std::endl;
237        }
238        if(newSolution != 0) delete [] newSolution;
239
240        int numberColumnCuts = theseCuts.sizeColCuts();
241        int numberRowCuts = theseCuts.sizeRowCuts();
242        violated = numberRowCuts + numberColumnCuts;
243       
244        //---------------------------------------------------------------------
245        // Apply column cuts
246        if (numberColumnCuts) {
247            double integerTolerance = getDblParam(AbcIntegerTolerance);
248            for (int i = 0; i < numberColumnCuts; ++i) {
249                const OsiColCut * thisCut = theseCuts.colCutPtr(i);
250                const CoinPackedVector & lbs = thisCut->lbs();
251                const CoinPackedVector & ubs = thisCut->ubs();
252                int j;
253                int n;
254                const int * which;
255                const double * values;
256                n = lbs.getNumElements();
257                which = lbs.getIndices();
258                values = lbs.getElements();
259                for (j = 0; j < n; ++j){
260                    int iColumn = which[j];
261                    double value = solution[iColumn];
262                    solver_->setColLower(iColumn, values[j]);
263                    if (value < values[j] - integerTolerance)
264                        violated = -1;   // violated, TODO: when happen?
265                    if (values[j] > upper[iColumn] + integerTolerance) {
266                        violated = -2;   // infeasible
267                        break;
268                    }
269                }
270                n = ubs.getNumElements();
271                which = ubs.getIndices();
272                values = ubs.getElements();
273                for (j = 0; j < n; ++j) {
274                    int iColumn = which[j];
275                    double value = solution[iColumn];
276                    solver_->setColUpper(iColumn, values[j]);
277                    if (value > values[j] + integerTolerance)
278                        violated = -1;
279                    if (values[j] < lower[iColumn] - integerTolerance) {
280                        violated = -2;   // infeasible
281                        break;
282                    }
283                }
284            }
285        }
286
287        if (violated == -2) {
288            feasible = false ;   
289            break ;    // break the cut generation loop
290        }
291
292        //---------------------------------------------------------------------
293        // Now apply the row (constraint) cuts.
294        int numberRowsNow = solver_->getNumRows();
295        assert(numberRowsNow == numberRowsAtStart + lastNumberCuts);
296        int numberToAdd = theseCuts.sizeRowCuts();
297        numberNewCuts = lastNumberCuts + numberToAdd;
298
299        // Get a basis by asking the solver for warm start information.
300        // Resize it (retaining the basis) so it can accommodate the cuts.
301        delete basis_;
302        basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart());
303        assert(basis_ != NULL); // make sure not volume
304        basis_->resize(numberRowsAtStart + numberNewCuts, numberColumns);
305
306        //  Now actually add the row cuts and reoptimise.
307        if (numberRowCuts > 0 || numberColumnCuts > 0) { 
308            if (numberToAdd > 0) { 
309                int i;
310                OsiRowCut * addCuts = new OsiRowCut [numberToAdd];
311                for (i = 0; i < numberToAdd; ++i) { 
312                    addCuts[i] = theseCuts.rowCut(i); 
313                }
314                solver_->applyRowCuts(numberToAdd, addCuts);
315                // AJK this caused a memory fault on Win32
316                delete [] addCuts;
317                for (i = 0; i < numberToAdd; ++i) { 
318                    cuts.insert(theseCuts.rowCut(i)); 
319                }
320                for (i = 0; i < numberToAdd; ++i) { 
321                    basis_->setArtifStatus(numberRowsNow + i,
322                                           CoinWarmStartBasis::basic); 
323                }
324                if (solver_->setWarmStart(basis_) == false) {
325                    throw CoinError("Fail setWarmStart() after cut install.",
326                                    "solveWithCuts", "SbbModel"); 
327                } 
328            }
329            feasible = resolve() ;
330        }
331        else { 
332            numberTries = 0; 
333        }
334
335        //---------------------------------------------------------------------
336        if (feasible) { 
337            int cutIterations = solver_->getIterationCount();
338            //takeOffCuts(cuts, whichGenerator,
339            // numberOldActiveCuts, numberNewCuts, true);
340            if (solver_->isDualObjectiveLimitReached()) { 
341                feasible = false;
342#ifdef ABC_DEBUG
343                double z = solver_->getObjValue();
344                double cut = getCutoff();
345                //      printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n",
346                //   z - cut, z, cut);
347#endif
348            }
349            if (feasible) { 
350                numberRowsAtStart = numberOldActiveCuts + 
351                    numberRowsAtContinuous_;
352                lastNumberCuts = numberNewCuts;
353                if ((direction * solver_->getObjValue() < 
354                    lastObjective + minimumDrop) &&  (numberPasses >= 3)) { 
355                    numberTries = 0; 
356                }
357                if (numberRowCuts+numberColumnCuts == 0 || cutIterations == 0)
358                { break; }
359                if (numberTries > 0) { 
360                    reducedCostFix();
361                    lastObjective = direction * solver_->getObjValue();
362                    lower = solver_->getColLower();
363                    upper = solver_->getColUpper();
364                    solution = solver_->getColSolution(); 
365                } 
366            } 
367        }
368
369        // We've lost feasibility
370        if (!feasible) { 
371            numberTries = 0;
372        }
373    } while (numberTries);
374    // END OF GENERATING CUTS
375
376    //------------------------------------------------------------------------
377    // Adjust the frequency of use for any of the cut generator
378    double thisObjective = solver_->getObjValue() * direction;
379    if (feasible && fullScan && numberCutGenerators_) {
380        double totalCuts = 0.0;
381        int i;
382        for (int i = 0; i < numberCutGenerators_; ++i) 
383        totalCuts += countRowCuts[i] + 5.0 * countColumnCuts[i];
384        // Root node or every so often - see what to turn off
385        if (!numberNodes_)
386            handler_->message(ABC_ROOT, messages_)
387                << numberNewCuts
388                << startObjective << thisObjective
389                << numberPasses
390                << CoinMessageEol;
391        int * count = new int[numberCutGenerators_];
392        memset(count, 0, numberCutGenerators_ * sizeof(int));
393        for (i = 0; i < numberNewCuts; ++i) 
394            count[whichGenerator[i]]++;
395        double small = (0.5 * totalCuts) / ((double) numberCutGenerators_);
396        for (i = 0; i < numberCutGenerators_; ++i) {
397            int howOften = generator_[i]->howOften();
398            if (howOften < -99)
399                continue;
400            if (howOften < 0 || howOften >= 1000000) {
401                // If small number switch mostly off
402                double thisCuts = countRowCuts[i] + 5.0 * countColumnCuts[i];
403                if (!thisCuts || howOften == -99) {
404                    if (howOften == -99)
405                        howOften = -100;
406                    else
407                        howOften = 1000000 + SCANCUTS; // wait until next time
408                } 
409                else if (thisCuts < small) {
410                    int k = (int) sqrt(small / thisCuts);
411                    howOften = k + 1000000;
412                } 
413                else {
414                    howOften = 1 + 1000000;
415                }
416            }
417            generator_[i]->setHowOften(howOften);
418            int newFrequency = generator_[i]->howOften() % 1000000;
419            // if (handler_->logLevel() > 1 || !numberNodes_)
420            if (!numberNodes_)
421                handler_->message(ABC_GENERATOR, messages_)
422                    << i
423                    << generator_[i]->cutGeneratorName()
424                    << countRowCuts[i]
425                    << countRowCuts[i] //<<count[i]
426                    << countColumnCuts[i]
427                    << newFrequency
428                    << CoinMessageEol;
429        }
430        delete [] count;
431    }
432
433    delete [] countRowCuts;
434    delete [] countColumnCuts;
435
436#ifdef CHECK_CUT_COUNTS
437    if (feasible) {
438        delete basis_;
439        basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart());
440        printf("solveWithCuts: Number of rows at end (only active cuts) %d\n",
441               numberRowsAtContinuous_+numberNewCuts+numberOldActiveCuts);
442        basis_->print(); 
443    }
444    if (numberNodes_ % 1000 == 0) {
445        messageHandler()->message(ABC_CUTS, messages_)
446            << numberNodes_
447            << numberNewCuts
448            << startObjective
449            << thisObjective
450            << numberPasses
451            << CoinMessageEol;
452    }
453#endif
454   
455
456    //takeOffCuts(cuts, whichGenerator, numberOldActiveCuts,
457    //          numberNewCuts, true);
458    incrementNodeCount();
459
460    return feasible;
461}
462
463//#############################################################################
464
465bool
466AbcModel::resolve()
467{
468    int iRow;
469    int numberRows = solver_->getNumRows();
470    const double * rowLower = solver_->getRowLower();
471    const double * rowUpper = solver_->getRowUpper();
472    bool feasible = true;
473    for (iRow = numberRowsAtContinuous_; iRow < numberRows; ++iRow) {
474        if (rowLower[iRow] > rowUpper[iRow] + 1.0e-8)
475            feasible = false;
476    }
477
478    // Reoptimize. Consider the possibility that we should fathom on bounds.
479    // But be careful --- where the objective takes on integral values, we may
480    // want to keep a solution where the objective is right on the cutoff.
481    if (feasible) { 
482        solver_->resolve();
483        numberIterations_ += getIterationCount();
484        feasible = (solver_->isProvenOptimal() &&
485                    !solver_->isDualObjectiveLimitReached()); 
486    }
487
488    return feasible;
489}
490
491//#############################################################################
492
493double 
494AbcModel::checkSolution (double cutoff, 
495                         const double *solution,
496                         bool fixVariables)
497{
498    return 0.0;
499}
500
501//#############################################################################
502
503bool
504AbcModel::setBestSolution(ABC_Message how,
505                          double & objectiveValue, 
506                          const double * solution, 
507                          bool fixVariables)
508{
509    double cutoff = getCutoff();
510    // Double check the solution to catch pretenders.
511    if (objectiveValue >= cutoff) {  // Bad news
512        if (objectiveValue > 1.0e30)
513            handler_->message(ABC_NOTFEAS1, messages_) << CoinMessageEol;
514        else
515            handler_->message(ABC_NOTFEAS2, messages_)
516                << objectiveValue << cutoff << CoinMessageEol;
517        return false;
518    }
519    else {  // Better solution
520        bestObjective_ = objectiveValue;
521        int numberColumns = solver_->getNumCols();
522        if (bestSolution_ == 0) {
523            bestSolution_ = new double[numberColumns];
524        }
525       
526        memcpy(bestSolution_, solution, numberColumns*sizeof(double));
527        cutoff = bestObjective_ - dblParam_[AbcCutoffIncrement];
528        setCutoff(cutoff);
529
530        if (how == ABC_ROUNDING)
531            numberHeuristicSolutions_++;
532        numberSolutions_++;
533//      std::cout << "cutoff = " << getCutoff()
534//                << "; objVal = " << bestObjective_
535//                << "; cutoffInc = " << dblParam_[AbcCutoffIncrement]
536//                << std::endl;
537       
538        handler_->message(how, messages_)
539            << bestObjective_ << numberIterations_
540            << numberNodes_
541            << CoinMessageEol;
542
543        return true;
544    }
545}
546
547//#############################################################################
548
549bool 
550AbcModel::feasibleSolution(int & numberIntegerInfeasibilities)
551{
552    bool feasible = true;
553    numberIntegerInfeasibilities = 0;
554    int i = -1;
555    const int numCols = getNumCols();
556
557    if (currentSolution_ != 0) {
558        delete [] currentSolution_;
559        currentSolution_ = 0;
560    }
561
562    currentSolution_ = new double [numCols];
563    memcpy(currentSolution_, solver_->getColSolution(),sizeof(double)*numCols);
564
565    for (i = 0; i < numberIntegers_; ++i) {
566        if ( ! checkInteger(currentSolution_[integerVariable_[i]]) ) {
567            ++numberIntegerInfeasibilities;
568            feasible = false;
569        }
570    }
571
572    return feasible;
573}
574
575//#############################################################################
576
577void 
578AbcModel::findIntegers(bool startAgain)
579{
580    assert(solver_);
581
582    int iColumn;   
583    int numberColumns = getNumCols();
584    const double *objCoeffs = getObjCoefficients();
585   
586    if (numberStrong_ == 0) {  // set up pseudocost list
587        pseudoList_ = new AbcPseudocost* [getNumCols()];
588        pseudoIndices_ = new int [getNumCols()];
589        for (iColumn = 0; iColumn < numberColumns; ++iColumn) {
590            pseudoList_[iColumn] = NULL;
591            pseudoIndices_[iColumn] = -1;
592        }
593    }
594
595    if (numberIntegers_ && !startAgain) return;
596    delete [] integerVariable_;
597    numberIntegers_ = 0;
598
599    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
600        if (isInteger(iColumn)) numberIntegers_++;
601    }
602
603    if (numberIntegers_) {
604        integerVariable_ = new int [numberIntegers_];
605        numberIntegers_=0;
606        for (iColumn = 0; iColumn < numberColumns; ++iColumn) {
607            if(isInteger(iColumn)) {
608                integerVariable_[numberIntegers_++] = iColumn;
609                if (numberStrong_ == 0) {
610                    double obj = fabs(objCoeffs[iColumn]);
611                    AbcPseudocost *pcost = new AbcPseudocost(iColumn,
612                                                             obj,
613                                                             0,
614                                                             obj,
615                                                             0);
616                    pseudoList_[iColumn] = pcost;
617                    //printf("numberIntegers_ = %d\n", numberIntegers_);
618                }
619                //printf("out\n");
620            }
621        }
622    } 
623    else {
624        handler_->message(ABC_NOINT, messages_) << CoinMessageEol ;
625    }
626
627   
628       
629}
630
631//#############################################################################
632// Add one generator
633void 
634AbcModel::addCutGenerator(CglCutGenerator * generator,
635                          int howOften, const char * name,
636                          bool normal, bool atSolution,
637                          bool whenInfeasible)
638{
639    AbcCutGenerator ** temp = generator_;
640    generator_ = new AbcCutGenerator * [numberCutGenerators_ + 1];
641    memcpy(generator_, temp, numberCutGenerators_*sizeof(AbcCutGenerator *));
642    delete[] temp;
643    generator_[numberCutGenerators_++]= 
644        new AbcCutGenerator(this, generator, howOften, name,
645                            normal, atSolution, whenInfeasible);
646                                                         
647}
648
649//#############################################################################
650// Add one heuristic
651void 
652AbcModel::addHeuristic(AbcHeuristic* generator)
653{
654  AbcHeuristic ** temp = heuristic_;
655  heuristic_ = new AbcHeuristic* [numberHeuristics_ + 1];
656  memcpy(heuristic_, temp, numberHeuristics_ * sizeof(AbcHeuristic *));
657  delete [] temp;
658  heuristic_[numberHeuristics_++] = generator;
659}
660
661//#############################################################################
662// Perform reduced cost fixing on integer variables. The variables in
663// question are already nonbasic at bound. We're just nailing down the
664// current situation.
665void AbcModel::reducedCostFix ()
666{ 
667    double cutoff = getCutoff();
668    double direction = solver_->getObjSense();
669    double gap = cutoff - solver_->getObjValue()*direction;
670    double integerTolerance = getDblParam(AbcIntegerTolerance);
671
672    const double* lower = solver_->getColLower();
673    const double* upper = solver_->getColUpper();
674    const double* solution = solver_->getColSolution();
675    const double* reducedCost = solver_->getReducedCost();
676
677    int numberFixed = 0 ;
678    for (int i = 0; i < numberIntegers_; i++) { 
679        int iColumn = integerVariable_[i];
680        double djValue = direction * reducedCost[iColumn];
681        if (upper[iColumn] - lower[iColumn] > integerTolerance) { 
682            if (solution[iColumn] < lower[iColumn] + integerTolerance && 
683                djValue > gap) { 
684                solver_->setColUpper(iColumn, lower[iColumn]);
685                numberFixed++; 
686            }
687            else if (solution[iColumn] > upper[iColumn] - integerTolerance && 
688                     -djValue > gap) { 
689                solver_->setColLower(iColumn, upper[iColumn]);
690                numberFixed++; 
691            } 
692        } 
693    } 
694}
695
696//#############################################################################
697#if 0
698void
699AbcModel::takeOffCuts(OsiCuts &newCuts, int *whichGenerator,
700                      int &numberOldActiveCuts, int &numberNewCuts,
701                      bool allowResolve)
702{
703    int firstOldCut = numberRowsAtContinuous_;
704    int totalNumberCuts = numberNewCuts + numberOldActiveCuts;
705    int *solverCutIndices = new int[totalNumberCuts];
706    int *newCutIndices = new int[numberNewCuts];
707    const CoinWarmStartBasis* ws;
708    CoinWarmStartBasis::Status status;
709    bool needPurge = true;
710   
711    // The outer loop allows repetition of purge in the event that
712    // reoptimisation changes the basis. To start an iteration, clear the
713    // deletion counts and grab the current basis.
714   
715    while (needPurge) {
716        int numberNewToDelete = 0;
717        int numberOldToDelete = 0;
718        int i;
719        ws = dynamic_cast<const CoinWarmStartBasis*>(solver_->getWarmStart());
720
721        // Scan the basis entries of the old cuts generated prior to this
722        // round of cut generation. Loose cuts are `removed'.
723        for (i = 0; i < numberOldActiveCuts; ++i) {
724            status = ws->getArtifStatus(i + firstOldCut);
725            if (status == CoinWarmStartBasis::basic) {
726                solverCutIndices[numberOldToDelete++] = i + firstOldCut;
727            }
728        }
729
730        // Scan the basis entries of the new cuts generated with this round
731        // of cut generation.  At this point, newCuts is the only record of
732        // the new cuts, so when we delete loose cuts from newCuts, they're
733        // really gone. newCuts is a vector, so it's most efficient to
734        // compress it (eraseRowCut) from back to front.
735        int firstNewCut = firstOldCut + numberOldActiveCuts;
736        int k = 0;
737        for (i = 0; i < numberNewCuts; ++i) {
738            status = ws->getArtifStatus(i + firstNewCut);
739            if (status == CoinWarmStartBasis::basic) {
740                solverCutIndices[numberNewToDelete + numberOldToDelete] =
741                    i + firstNewCut ;
742                newCutIndices[numberNewToDelete++] = i;
743            }
744            else { // save which generator did it
745                whichGenerator[k++] = whichGenerator[i];
746            }
747        }
748        for (i = numberNewToDelete - 1 ; i >= 0 ; i--) {
749            int iCut = newCutIndices[i];
750            newCuts.eraseRowCut(iCut);
751        }
752
753        // Did we delete anything? If so, delete the cuts from the constraint
754        // system held in the solver and reoptimise unless we're forbidden
755        // to do so. If the call to resolve() results in pivots, there's the
756        // possibility we again have basic slacks. Repeat the purging loop.
757
758        if (numberNewToDelete  + numberOldToDelete > 0) {
759            solver_->deleteRows(numberNewToDelete + numberOldToDelete,
760                                solverCutIndices);
761            numberNewCuts -= numberNewToDelete;
762            numberOldActiveCuts -= numberOldToDelete;
763#           ifdef ABC_DEBUG
764            std::cout << "takeOffCuts: purged " << numberOldToDelete << "+"
765                      << numberNewToDelete << " cuts." << std::endl;
766#           endif
767            if (allowResolve) {
768                solver_->resolve();
769                if (solver_->getIterationCount() == 0) {
770                    needPurge = false;
771                }
772#           ifdef ABC_DEBUG
773                else {
774                    std::cout << "Repeating purging loop. "
775                              << solver_->getIterationCount() << " iters."
776                              << std::endl;
777                }
778#           endif
779            }
780            else {
781                needPurge = false;
782            }
783        }
784        else {
785            needPurge = false;
786        }
787    }
788
789    delete ws;
790    delete [] solverCutIndices;
791    delete [] newCutIndices;
792}
793#endif
794
795//#############################################################################
796#if 1
797//void
798//AbcModel::takeOffCuts(OsiCuts &newCuts, int *whichGenerator,
799//                    int &numberOldActiveCuts, int &numberNewCuts,
800//                    bool allowResolve)
801void
802AbcModel::takeOffCuts()
803{
804//    assert(!numberOldActiveCuts);
805    int totalNumberCuts = solver()->getNumRows() - numberRowsAtContinuous_;
806    int *solverCutIndices = new int[totalNumberCuts];
807    //  const CoinWarmStartBasis* ws;
808   
809    for (int i = 0; i < totalNumberCuts; ++i) {
810        solverCutIndices[i] = i + numberRowsAtContinuous_;
811    }
812   
813    // Delete all new cuts
814    solver_->deleteRows(totalNumberCuts, solverCutIndices);
815    solver_->setWarmStart(sharedBasis_);
816    //numberOldActiveCuts = numberNewCuts = 0;
817    delete []  solverCutIndices;
818}
819#endif
820
821//#############################################################################
822/**
823  This routine sets the objective cutoff value used for fathoming and
824  determining monotonic variables.
825
826  If the fathoming discipline is strict, a small tolerance is added to the
827  new cutoff. This avoids problems due to roundoff when the target value
828  is exact. The common example would be an IP with only integer variables in
829  the objective. If the target is set to the exact value z of the optimum,
830  it's possible to end up fathoming an ancestor of the solution because the
831  solver returns z+epsilon.
832
833  Determining if strict fathoming is needed is best done by analysis.
834  In sbb, that's analyseObjective. The default is false.
835
836  In sbb we always minimize so add epsilon
837*/
838void AbcModel::setCutoff (double value)
839{ 
840    double tol = 0;
841    int fathomStrict = getIntParam(AbcFathomDiscipline);
842    double direction = solver_->getObjSense();
843    if (fathomStrict == 1) { 
844        solver_->getDblParam(OsiDualTolerance, tol);
845        tol = tol * (1 + fabs(value));
846        value += tol; 
847    }
848   
849    // Solvers know about direction
850    solver_->setDblParam(OsiDualObjectiveLimit, value * direction); 
851}
852
853//#############################################################################
854// Initial solve and find integers
855bool
856AbcModel::setupSelf()
857{   
858    bool feasible = true;
859    solver_->messageHandler()->setLogLevel(0);
860    initialSolve();
861    sharedBasis_ = dynamic_cast<CoinWarmStartBasis*>
862        (solver_->getWarmStart());
863   
864# ifdef ABC_DEBUG_MORE
865    std::string problemName;
866    solver_->getStrParam(OsiProbName, problemName);
867    printf("Problem name - %s\n", problemName.c_str());
868    solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0);
869# endif
870
871    status_ = 0;
872
873    findIntegers(true);
874
875    bestObjective_ = 1.0e50;
876    double direction = solver_->getObjSense();
877
878    int numberColumns = getNumCols();
879    if (!currentSolution_)
880        currentSolution_ = new double[numberColumns];
881
882    continuousSolver_ = solver_->clone();
883    numberRowsAtContinuous_ = getNumRows();
884
885    maximumNumberCuts_ = 0;
886    currentNumberCuts_ = 0;
887
888    // FIXME:
889   
890    return feasible;
891}
892
893//#############################################################################
894// Send model and root so that initial solve
895AlpsEncoded* 
896AbcModel::encode() const 
897{ 
898    AlpsReturnCode status = ALPS_OK;
899
900    AlpsEncoded* encoded = new AlpsEncoded("ALPS_MODEL");
901
902    //------------------------------------------------------
903    // Encode Alps part.
904    //------------------------------------------------------
905
906    status = encodeAlps(encoded);
907   
908    //------------------------------------------------------
909    // Encode Abc part.
910    //------------------------------------------------------
911
912    // Write the model data into representation_
913    const CoinPackedMatrix* matrixByCol = solver_->getMatrixByCol();
914    int numRows = getNumRows();
915    encoded->writeRep(numRows);
916    int numCols = getNumCols();
917    encoded->writeRep(numCols);
918#if defined(ABC_DEBUG_MORE)
919    std::cout << "AbcModel::encode()-- numRows="<< numRows << "; numCols=" 
920              << numCols << std::endl;
921#endif
922
923    const double* collb = solver_->getColLower();
924    encoded->writeRep(collb, numCols);
925    const double* colub = solver_->getColUpper();
926    encoded->writeRep(colub, numCols);
927    const double* obj = solver_->getObjCoefficients();
928    encoded->writeRep(obj, numCols);
929    const double objSense = solver_->getObjSense();
930    encoded->writeRep(objSense);
931    const double* rowlb = solver_->getRowLower();
932    encoded->writeRep(rowlb, numRows);
933    const double* rowub = solver_->getRowUpper();
934    encoded->writeRep(rowub, numRows);
935    int numElements = solver_->getNumElements();
936    encoded->writeRep(numElements);
937    const double* elementValue = matrixByCol->getElements();
938    encoded->writeRep(elementValue, numElements);
939    const CoinBigIndex* colStart = matrixByCol->getVectorStarts();
940    int numStart = numCols + 1;
941    encoded->writeRep(colStart, numStart);
942    const int* index = matrixByCol->getIndices();
943    encoded->writeRep(index, numElements);
944    encoded->writeRep(numberIntegers_);
945    encoded->writeRep(integerVariable_, numberIntegers_);
946#if defined(ABC_DEBUG_MORE)
947    std::cout << "AbcModel::encode()-- objSense="<< objSense
948              << "; numElements="<< numElements
949              << "; numberIntegers_=" << numberIntegers_
950              << "; numStart = " << numStart <<std::endl;
951#endif
952#if defined(ABC_DEBUG_MORE)
953    std::cout << "rowub=";
954    for (int i = 0; i < numRows; ++i){
955        std::cout <<rowub[i]<<" ";
956    }
957    std::cout << std::endl;
958    std::cout << "elementValue=";
959    for (int j = 0; j < numElements; ++j) {
960        std::cout << elementValue[j] << " ";
961    }
962    std::cout << std::endl;   
963#endif
964
965    return encoded;
966}
967
968//#############################################################################
969// Decode and load model data to LP solver.
970void
971AbcModel::decodeToSelf(AlpsEncoded& encoded) 
972{
973    AlpsReturnCode status = ALPS_OK;
974
975    //------------------------------------------------------
976    // Decode Alps part.
977    //------------------------------------------------------
978
979    status = decodeAlps(encoded);
980
981    //------------------------------------------------------
982    // Decode Abc part.
983    //------------------------------------------------------
984
985    int numRows;
986    encoded.readRep(numRows);
987    int numCols;
988    encoded.readRep(numCols);   
989#if defined(ABC_DEBUG_MORE)
990    std::cout << "AbcModel::decode()-- numRows="<< numRows << "; numCols=" 
991              << numCols << std::endl;
992#endif
993    double* collb;
994    encoded.readRep(collb, numCols);
995    double* colub;
996    encoded.readRep(colub, numCols);
997    double* obj;
998    encoded.readRep(obj, numCols);
999    double objSense;
1000    encoded.readRep(objSense);
1001    double* rowlb;
1002    encoded.readRep(rowlb, numRows);
1003    double* rowub;
1004    encoded.readRep(rowub, numRows);
1005    int numElements;
1006    encoded.readRep(numElements);
1007    double* elementValue;
1008    encoded.readRep(elementValue, numElements);
1009    CoinBigIndex* colStart;
1010    int numStart = numCols + 1;
1011    encoded.readRep(colStart, numStart);
1012    int* index;
1013    encoded.readRep(index, numElements);
1014    encoded.readRep(numberIntegers_);
1015    encoded.readRep(integerVariable_, numberIntegers_);
1016#if defined(ABC_DEBUG_MORE)
1017    std::cout << "AbcModel::decode()-- objSense="<< objSense
1018              <<  "; numElements="<< numElements
1019              << "; numberIntegers_=" << numberIntegers_
1020              << "; numStart = " << numStart <<std::endl;
1021#endif
1022#if defined(ABC_DEBUG_MORE)
1023    std::cout << "rowub=";
1024    for (int i = 0; i < numRows; ++i){
1025        std::cout <<rowub[i]<<" ";
1026    }
1027    std::cout << std::endl;
1028    std::cout << "elementValue=";
1029    for (int j = 0; j < numElements; ++j) {
1030        std::cout << elementValue[j] << " ";
1031    }
1032    std::cout << std::endl; 
1033    std::cout << "index=";
1034    for (int j = 0; j < numElements; ++j) {
1035        std::cout << index[j] << " ";
1036    }
1037    std::cout << std::endl; 
1038    std::cout << "colStart=";
1039    for (int j = 0; j < numElements+1; ++j) {
1040        std::cout << colStart[j] << " ";
1041    }
1042    std::cout << std::endl;   
1043#endif
1044
1045    // Check if solver_ is declared in main
1046    assert(solver_);
1047
1048    //-------------------------------------------------------------------------
1049    // load the standardized problem into stdSi   
1050#if 0  // load matrix doesn't work. Don't know why.
1051    CoinPackedMatrix * matrixByCol =
1052        new CoinPackedMatrix(true, numCols, numRows, numElements,
1053                             elementValue, index, colStart, 0);
1054    solver_->loadProblem(*matrixByCol,
1055                         collb, colub,   
1056                         obj,
1057                         rowlb, rowub);
1058#endif
1059    //-------------------------------------------------------------------------
1060    solver_->loadProblem(numCols, numRows,
1061                         colStart, index, elementValue,
1062                         collb, colub, 
1063                         obj,
1064                         rowlb, rowub);
1065
1066    solver_->setObjSense(objSense);
1067    solver_->setInteger(integerVariable_, numberIntegers_);
1068
1069    delete [] collb;
1070    collb = NULL;
1071    delete [] colub;
1072    colub = NULL;
1073    delete [] obj;
1074    obj = NULL;
1075    delete [] rowlb;
1076    rowlb = NULL;
1077    delete [] rowub;
1078    rowub = NULL;
1079    delete [] elementValue;
1080    elementValue = NULL;
1081    delete [] colStart;
1082    colStart = NULL;
1083    delete [] index;
1084    index = NULL;
1085}
Note: See TracBrowser for help on using the repository browser.