source: trunk/Alps/examples/Abc/AbcHeuristic.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: 12.3 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 SbbHeuristic.cpp
20//#############################################################################
21
22#include <cassert>
23#include <cmath>
24#include <cfloat>
25
26#include "OsiSolverInterface.hpp"
27#include "AbcModel.h"
28//#include "AbcMessage.hpp"
29#include "AbcHeuristic.h"
30
31// Default Constructor
32AbcHeuristic::AbcHeuristic() 
33    :model_(NULL)
34{
35}
36
37// Constructor from model
38AbcHeuristic::AbcHeuristic(AbcModel & model)
39{
40    model_ = &model; 
41}
42
43// Destructor
44AbcHeuristic::~AbcHeuristic ()
45{
46}
47
48// update model
49void AbcHeuristic::setModel(AbcModel * model)
50{
51    model_ = model;
52}
53
54// Default Constructor
55AbcRounding::AbcRounding() 
56    :AbcHeuristic()
57{
58    // matrix and row copy will automatically be empty
59}
60
61// Constructor from model
62AbcRounding::AbcRounding(AbcModel & model)
63    :AbcHeuristic(model)
64{
65    assert(model.solver());
66}
67
68// Destructor
69AbcRounding::~AbcRounding ()
70{
71}
72
73// Clone
74AbcHeuristic *
75AbcRounding::clone() const
76{
77    return new AbcRounding(*this);
78}
79
80// Copy constructor
81AbcRounding::AbcRounding(const AbcRounding & rhs)
82    :
83    AbcHeuristic(rhs),
84    matrix_(rhs.matrix_),
85    matrixByRow_(rhs.matrixByRow_),
86    seed_(rhs.seed_)
87{
88}
89// See if rounding will give solution
90// Sets value of solution
91// Assumes rhs for original matrix still okay
92// At present only works with integers
93// Fix values if asked for
94// Returns 1 if solution, 0 if not
95int
96AbcRounding::solution(double & solutionValue,
97                      double * betterSolution)
98{
99    // Get a copy of original matrix (and by row for rounding);
100    matrix_ = *(model_->solver()->getMatrixByCol());
101    matrixByRow_ = *(model_->solver()->getMatrixByRow());
102    seed_=1;
103
104    OsiSolverInterface * solver = model_->solver();
105    const double * lower = solver->getColLower();
106    const double * upper = solver->getColUpper();
107    const double * rowLower = solver->getRowLower();
108    const double * rowUpper = solver->getRowUpper();
109    const double * solution = solver->getColSolution();
110    const double * objective = solver->getObjCoefficients();
111    double integerTolerance = 1.0e-5;
112    //model_->getDblParam(AbcModel::AbcIntegerTolerance);
113    double primalTolerance;
114    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
115
116    int numberRows = matrix_.getNumRows();
117
118    int numberIntegers = model_->numberIntegers();
119    const int * integerVariable = model_->integerVariable();
120    int i;
121    double direction = solver->getObjSense();
122    double newSolutionValue = direction * solver->getObjValue();
123    int returnCode = 0;
124
125    // Column copy
126    const double * element = matrix_.getElements();
127    const int * row = matrix_.getIndices();
128    const int * columnStart = matrix_.getVectorStarts();
129    const int * columnLength = matrix_.getVectorLengths();
130    // Row copy
131    const double * elementByRow = matrixByRow_.getElements();
132    const int * column = matrixByRow_.getIndices();
133    const int * rowStart = matrixByRow_.getVectorStarts();
134    const int * rowLength = matrixByRow_.getVectorLengths();
135
136    // Get solution array for heuristic solution
137    int numberColumns = solver->getNumCols();
138    double * newSolution = new double [numberColumns];
139    memcpy(newSolution, solution, numberColumns * sizeof(double));
140
141    double * rowActivity = new double[numberRows];
142    memset(rowActivity, 0, numberRows*sizeof(double));
143    for (i = 0; i < numberColumns; i++) {
144        int j;
145        double value = newSolution[i];
146        if (value) {
147            for (j = columnStart[i];
148                 j < columnStart[i] + columnLength[i]; j++) {
149                int iRow = row[j];
150                rowActivity[iRow] += value*element[j];
151            }
152        }
153    }
154    // check was feasible - if not adjust (cleaning may move)
155    for (i = 0; i < numberRows; i++) {
156        if(rowActivity[i] < rowLower[i]) {
157            //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
158            rowActivity[i] = rowLower[i];
159        } else if(rowActivity[i] > rowUpper[i]) {
160            //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
161            rowActivity[i] = rowUpper[i];
162        }
163    }
164    for (i = 0; i < numberIntegers; i++) {
165        int iColumn = integerVariable[i];
166        double value = newSolution[iColumn];
167        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
168            double below = floor(value);
169            double newValue = newSolution[iColumn];
170            double cost = direction * objective[iColumn];
171            double move;
172            if (cost > 0.0) {
173                // try up
174                move = 1.0 - (value - below);
175            } else if (cost < 0.0) {
176                // try down
177                move = below - value;
178            } else {
179                // won't be able to move unless we can grab another variable
180                // just for now go down
181                move = below-value;
182            }
183            newValue += move;
184            newSolution[iColumn] = newValue;
185            newSolutionValue += move * cost;
186            int j;
187            for (j = columnStart[iColumn];
188                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
189                int iRow = row[j];
190                rowActivity[iRow] += move * element[j];
191            }
192        }
193    }
194
195    double penalty = 0.0;
196 
197    // see if feasible
198    for (i = 0; i < numberRows; i++) {
199        double value = rowActivity[i];
200        double thisInfeasibility = 0.0;
201        if (value < rowLower[i] - primalTolerance)
202            thisInfeasibility = value - rowLower[i];
203        else if (value > rowUpper[i] + primalTolerance)
204            thisInfeasibility = value - rowUpper[i];
205        if (thisInfeasibility) {
206            // See if there are any slacks I can use to fix up
207            // maybe put in coding for multiple slacks?
208            double bestCost = 1.0e50;
209            int k;
210            int iBest = -1;
211            double addCost = 0.0;
212            double newValue = 0.0;
213            double changeRowActivity = 0.0;
214            double absInfeasibility = fabs(thisInfeasibility);
215            for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
216                int iColumn = column[k];
217                if (columnLength[iColumn] == 1) {
218                    double currentValue = newSolution[iColumn];
219                    double elementValue = elementByRow[k];
220                    double lowerValue = lower[iColumn];
221                    double upperValue = upper[iColumn];
222                    double gap = rowUpper[i] - rowLower[i];
223                    double absElement = fabs(elementValue);
224                    if (thisInfeasibility * elementValue > 0.0) {
225                        // we want to reduce
226                        if ((currentValue - lowerValue) * absElement >= 
227                            absInfeasibility) {
228                           
229                            // possible - check if integer
230                            double distance = absInfeasibility / absElement;
231                            double thisCost = 
232                                -direction * objective[iColumn] * distance;
233                            if (solver->isInteger(iColumn)) {
234                                distance = ceil(distance - primalTolerance);
235                                assert (currentValue - distance >= 
236                                        lowerValue - primalTolerance);
237                                if (absInfeasibility - distance * absElement
238                                    < -gap - primalTolerance)
239                                    thisCost = 1.0e100; // no good
240                                else
241                                    thisCost = 
242                                        -direction*objective[iColumn]*distance;
243                            }
244                            if (thisCost < bestCost) {
245                                bestCost = thisCost;
246                                iBest = iColumn;
247                                addCost = thisCost;
248                                newValue = currentValue - distance;
249                                changeRowActivity = -distance * elementValue;
250                            }
251                        }
252                    } else {
253                        // we want to increase
254                        if ((upperValue - currentValue) * absElement >= 
255                            absInfeasibility) {
256                            // possible - check if integer
257                            double distance = absInfeasibility / absElement;
258                            double thisCost = 
259                                direction * objective[iColumn] * distance;
260                            if (solver->isInteger(iColumn)) {
261                                distance = ceil(distance - 1.0e-7);
262                                assert (currentValue - distance <= 
263                                        upperValue + primalTolerance);
264                                if (absInfeasibility - distance * absElement
265                                    < -gap - primalTolerance)
266                                    thisCost = 1.0e100; // no good
267                                else
268                                    thisCost = 
269                                        direction*objective[iColumn]*distance;
270                            }
271                            if (thisCost < bestCost) {
272                                bestCost = thisCost;
273                                iBest = iColumn;
274                                addCost = thisCost;
275                                newValue = currentValue + distance;
276                                changeRowActivity = distance * elementValue;
277                            }
278                        }
279                    }
280                }
281            }
282            if (iBest >= 0) {
283                /*printf("Infeasibility of %g on row %d cost %g\n",
284                  thisInfeasibility,i,addCost);*/
285                newSolution[iBest] = newValue;
286                thisInfeasibility = 0.0;
287                newSolutionValue += addCost;
288                rowActivity[i] += changeRowActivity;
289            }
290            penalty += fabs(thisInfeasibility);
291        }
292    }
293
294    // Could also set SOS (using random) and repeat
295    if (!penalty) {
296        // See if we can do better
297        //seed_++;
298        //CoinSeedRandom(seed_);
299        // Random number between 0 and 1.
300        double randomNumber = CoinDrand48();
301        int iPass;
302        int start[2];
303        int end[2];
304        int iRandom = (int) (randomNumber * ((double) numberIntegers));
305        start[0] = iRandom;
306        end[0] = numberIntegers;
307        start[1] = 0;
308        end[1] = iRandom;
309        for (iPass = 0; iPass < 2; iPass++) {
310            int i;
311            for (i = start[iPass]; i < end[iPass]; i++) {
312                int iColumn = integerVariable[i];
313                double value = newSolution[iColumn];
314                assert(fabs(floor(value + 0.5) - value) < integerTolerance);
315                double cost = direction * objective[iColumn];
316                double move = 0.0;
317                if (cost > 0.0)
318                    move = -1.0;
319                else if (cost < 0.0)
320                    move = 1.0;
321                while (move) {
322                    bool good = true;
323                    double newValue = newSolution[iColumn] + move;
324                    if (newValue < lower[iColumn] - primalTolerance||
325                        newValue > upper[iColumn] + primalTolerance) {
326                        move = 0.0;
327                    } else {
328                        // see if we can move
329                        int j;
330                        for (j = columnStart[iColumn];
331                             j < columnStart[iColumn] + columnLength[iColumn];
332                             j++) {
333                            int iRow = row[j];
334                            double newActivity = 
335                                rowActivity[iRow] + move*element[j];
336                            if (newActivity < rowLower[iRow] - primalTolerance
337                                ||
338                                newActivity > rowUpper[iRow]+primalTolerance) {
339                                good = false;
340                                break;
341                            }
342                        }
343                        if (good) {
344                            newSolution[iColumn] = newValue;
345                            newSolutionValue += move * cost;
346                            int j;
347                            for (j = columnStart[iColumn];
348                                 j < columnStart[iColumn] + 
349                                     columnLength[iColumn]; j++) {
350                                int iRow = row[j];
351                                rowActivity[iRow] += move*element[j];
352                            }
353                        } else {
354                            move=0.0;
355                        }
356                    }
357                }
358            }
359        }
360        if (newSolutionValue < solutionValue) {
361            // paranoid check
362            memset(rowActivity, 0, numberRows * sizeof(double));
363            for (i = 0; i < numberColumns; i++) {
364                int j;
365                double value = newSolution[i];
366                if (value) {
367                    for (j = columnStart[i];
368                         j < columnStart[i] + columnLength[i]; j++) {
369                        int iRow = row[j];
370                        rowActivity[iRow] += value * element[j];
371                    }
372                }
373            }
374            // check was approximately feasible
375            bool feasible = true;
376            for (i = 0; i < numberRows; i++) {
377                if(rowActivity[i] < rowLower[i]) {
378                    if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
379                        feasible = false;
380                } else if(rowActivity[i] > rowUpper[i]) {
381                    if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
382                        feasible = false;
383                }
384            }
385            if (feasible) {
386                // new solution
387                memcpy(betterSolution, newSolution, 
388                       numberColumns * sizeof(double));
389                solutionValue = newSolutionValue;
390                //printf("** Solution of %g found by rounding\n",newSolutionValue);
391                returnCode=1;
392            } else {
393                // Can easily happen
394                //printf("Debug AbcRounding giving bad solution\n");
395            }
396        }
397    }
398    delete [] newSolution;
399    delete [] rowActivity;
400    return returnCode;
401}
402// update model
403void AbcRounding::setModel(AbcModel * model)
404{
405    model_ = model;
406    // Get a copy of original matrix (and by row for rounding);
407    assert(model_->solver());
408    matrix_ = *model_->solver()->getMatrixByCol();
409    matrixByRow_ = *model_->solver()->getMatrixByRow();
410}
Note: See TracBrowser for help on using the repository browser.