source: pico/trunk/src/milp/pico/feasibilityPump.cpp @ 5603

Revision 5603, 31.6 KB checked in by jwatson, 6 years ago (diff)

Numerous changes/improvements related to cycle detection and escape from local minima.
Configured the current EH variant to use triangular merit functions, to verify FP-equivalence.
When verified, I will be making the merit function type a command-line argument.

Line 
1/*  _________________________________________________________________________
2 *
3 *  PICO: A C++ library of scalable branch-and-bound and related methods.
4 *  Copyright (c) 2007, Sandia National Laboratories.
5 *  This software is distributed under the GNU Lesser General Public License.
6 *  For more information, see the README file in the top PICO directory.
7 *  _________________________________________________________________________
8 */
9
10
11#include <acro_config.h>
12#include <pebbl/gRandom.h>
13#include <pico/milp.h>
14#include <pico/milpNode.h>
15
16#include <pico/feasibilityPump.h>
17
18#include <CoinPackedVector.hpp>
19#include <CoinWarmStart.hpp>
20
21#include <iostream>
22#include <iomanip>
23#include <limits>
24
25#include <math.h>
26#include <assert.h>
27
28// comparison function used in std::sort call found below.
29// sorts in decreasing order.
30static bool compareDists(const std::pair<double,int> & a,
31                         const std::pair<double,int> & b)
32{
33  return (a.first > b.first);
34}
35
36// a pretty-print / debug utility.
37template <typename X>
38void prettyPrint(const char * prefix,
39                 const X * vector_data,
40                 unsigned int vector_length,
41                 std::ostream &ostr,
42                 int width = 4)
43{
44  ostr << prefix << "=" << std::endl;
45  ostr << "[";
46  for (unsigned int i = 0; i < vector_length; ++i)
47    {
48      ostr.precision(2);
49      ostr << std::setw(width) << vector_data[i];
50      if (i != vector_length - 1)
51        {
52          ostr << " ";
53        }
54    }
55  ostr << "]" << std::endl;
56}
57
58namespace pico
59{
60
61  feasibilityPump::feasibilityPump(void)
62  {
63    name = "Feasibility pump";
64
65    FPDebug = 0;
66    ParameterSet::create_categorized_parameter("FPDebug",FPDebug,
67                                               "<int>","0",
68                                               "Specific debug level for feasibility pump",
69                                               "Debugging",
70                                               utilib::ParameterPositive<int>());
71
72    FPMetaIterations = 1;
73    ParameterSet::create_categorized_parameter("FPMetaIterations",FPMetaIterations,
74                                               "<int>","1",
75                                               "Feasibility pump meta-iteration limit",
76                                               "Feasibility Pump Heuristic",
77                                               utilib::ParameterPositive<int>());
78
79    FPRootIterations = 2000;
80    ParameterSet::create_categorized_parameter("FPRootIterations",FPRootIterations,
81                                               "<int>","2000",
82                                               "Limit on the cumulative feasibility pump iterations at the root node.",
83                                               "Feasibility Pump Heuristic",
84                                               utilib::ParameterPositive<int>());
85
86    FPInteriorIterations = 50;
87    ParameterSet::create_categorized_parameter("FPInteriorIterations",FPInteriorIterations,
88                                               "<int>","50",
89                                               "Limit on the cumulative feasibility pump iterations at interior nodes.",
90                                               "Feasibility Pump Heuristic",
91                                               utilib::ParameterPositive<int>());
92
93    FPRunTime = 300;
94    ParameterSet::create_categorized_parameter("FPRunTime",FPRunTime,
95                                               "<double>","300",
96                                               "Limit on the cumulative CPU seconds spent in the feasibility pump heuristic",
97                                               "Feasibility Pump Heuristic",
98                                               utilib::ParameterPositive<double>());
99
100    FPRandomInitialSolution = false;
101    ParameterSet::create_categorized_parameter("FPRandomInitialSolution",FPRandomInitialSolution,
102                                               "<bool>","false",
103                                               "Starts feasibilty pump from a random solution",
104                                               "Feasibility Pump Heuristic");
105
106    FPImproveIncumbent = true;
107    ParameterSet::create_categorized_parameter("FPImproveIncumbent",FPImproveIncumbent,
108                                               "<bool>","true",
109                                               "Allows feasibility pump to continue inner-loop iterations in attempt to improve an initial incumbent",
110                                               "Feasibility Pump Heuristic");
111
112    FPImprovementAlpha = 0.50;
113    ParameterSet::create_categorized_parameter("FPImprovementAlpha",FPImprovementAlpha,
114                                               "<double>","0.50",
115                                               "When improving incumbents, the fraction of the absolute gap used to form the new objective cut",
116                                               "Feasibility Pump Heuristic");
117
118    FPMeanFlipCount = 20;
119    ParameterSet::create_categorized_parameter("FPMeanFlipCount",FPMeanFlipCount,
120                                               "<int>","20",
121                                               "Mean number of variables 'flipped' per feasibility pump iteration, if rounded solution remains unchanged",
122                                               "Feasibility Pump Heuristic",
123                                               utilib::ParameterPositive<int>());
124
125    FPRoundingType = 0;
126    ParameterSet::create_categorized_parameter("FPRoundingType",FPRoundingType,
127                                               "<int>","0",
128                                               "The type of rounding used by FP to form a new x~ target vector from the current LP solution x*. 0 indicates standard deterministic. 1 indicates randomized rounding around the interval [0.33, 0.66]",
129                                               "Feasibility Pump Heuristic",
130                                               utilib::ParameterPositive<int>());
131
132    FPPerturbationInterval = 100;
133    ParameterSet::create_categorized_parameter("FPPerturbationInterval",FPPerturbationInterval,
134                                               "<int>","100",
135                                               "The number of FP iterations between random perturbations of the current x~ solution",
136                                               "Feasibility Pump Heuristic",
137                                               utilib::ParameterPositive<int>());
138
139  }
140
141  feasibilityPump::~feasibilityPump(void)
142  {
143    // no dynamic memory to worry about.
144  }
145
146  void feasibilityPump::reset(MILP* mGlobal_, lpShareObj* lpShare_)
147  {
148    lpBasedHeuristic::reset(mGlobal_, lpShare_);
149
150    // Serial parameter settings that are hierarchical/not automatic
151
152    if (parameter_initialized("FPDebug"))
153      {
154        setDebug(FPDebug);
155      }
156
157    // FP needs the objective cut when dealing with
158    // non-root-node solves.
159    useObjCut = true;
160  }
161
162  void feasibilityPump::setupLP(MILPNode * sp)
163  {
164    // we currently execute FP whenever the heuristic management framework tells us to...
165    currentDepth = sp->depth;
166    executeFP = true;
167
168    DEBUGPR(10, ucout << "Setting up feasibility pump to execute at depth=" << sp->depth << std::endl << Flush);
169
170    lpBasedHeuristic::setupLP(sp);
171  }
172
173  MILPSolution * feasibilityPump::runHeuristic(void)
174  {
175    // execute FP according to whatever rule is encoded in the setupLP method.
176    if (executeFP == false)
177      {
178        return 0;
179      }
180
181    // stops the LP from thinking that it's just solving the relaxation -
182    // we  need to know which variables are integer.
183    // IMPT: have to un-do on the way out.
184    lp->setCutting();
185
186    // verify that the rounding type is set to a sane value.
187    if (FPRoundingType >= 2)
188      {
189        DEBUGPR(10, ucout << "***WARNING: FPRoundingType parameter has an unknown value=" << FPRoundingType << std::endl << Flush);
190        return 0;
191      }
192
193    DEBUGPR(100, ucout << "-----------------------------------------------" << std::endl << Flush);
194    DEBUGPR(10, ucout << "Executing feasibility pump" << std::endl << Flush);
195
196    double start_time = CPUSeconds();
197
198    int cumulative_iterations(0);
199    int allowable_iterations;
200    if (currentDepth == 1) // depth is 1-based in PICO/PEBBL.
201      {
202        allowable_iterations = FPRootIterations;
203      }
204    else
205      {
206        allowable_iterations = FPInteriorIterations;
207      }
208
209    int num_vars(lp->getNumCols());
210
211    // verify that any integer variables are in fact binary - the current
212    // code doesn't deal with the general integer case, as that's more
213    // than a bit of pain.
214    for (int i = 0; i < num_vars; ++i)
215      {
216        if ((lp->isContinuous(i) == false) &&
217            (lp->isBinary(i) == false))
218          {
219            DEBUGPR(10, ucout << "Not all integer variables are binary - "
220                    "feasibility pump is currently only coded for binary MIPs"
221                    << std::endl << Flush);
222            return 0;
223          }
224      }
225
226    DEBUGPR(50, ucout << "Limit on number of feasibility pump meta iterations="
227            << FPMetaIterations << std::endl << Flush);
228    DEBUGPR(50, ucout << "Limit on cumulative number of feasibility pump "
229            "iterations=" << allowable_iterations << std::endl << Flush);
230    DEBUGPR(50, ucout << "Limit (in seconds) on the cumulative feasibility pump "
231            " cpu time=" << FPRunTime << std::endl << Flush);
232
233    // statistics for the best-so-far incumbent.
234    double best_objective;
235    DoubleVector best_solution(num_vars);
236    if (lp->getObjSense() == pebbl::minimize)
237      {
238        best_objective = std::numeric_limits<double>::max();
239      }
240    else
241      {
242        best_objective = std::numeric_limits<double>::min();
243      }
244
245    int executed_iterations(0);
246
247    for (int meta_iteration = 1;
248         meta_iteration <= FPMetaIterations;
249         ++meta_iteration)
250      {
251        DEBUGPR(10, ucout << "Executing meta-level iteration="
252                << meta_iteration << std::endl << Flush);
253
254        ++executed_iterations;
255
256        // clone the LP interface to solve the feasibility pump sub-problems.
257        // NOTE: if you do this in the outer loop, i.e., create it once,
258        //       then re-solves with the LP* relaxation initial solution
259        //       immediately returned with a feasible (high-quality) solution.
260        //       this wasn't happening for random initial solutions. not
261        //       resolved yet as to why this happens...
262        // NOTE: I moved the clone/resolve to this point, because the
263        //       lp object has no objective/solution - for reasons we
264        //       don't understand at all.
265        OsiSolverInterface * fp_solver_interface = lp->OsiOnlyCopy();
266        // Cbc does this - isn't clear to me why, but...
267        fp_solver_interface->resolve();
268
269        // depending on the instance and the time allocated to FP, the
270        // resolve may take more time that is allocated to FP. in this
271        // case, break out of the meta-iteration loop and terminate.
272        if ((CPUSeconds() - start_time) >= FPRunTime)
273          {
274            DEBUGPR(10, ucout << "Time limit exceeded during resolve() call to cloned LP solver for FP" << std::endl << Flush);
275            break;
276          }
277
278        // form an integer (but likely infeasible) starting solution.
279        double * initial_solution = new double[num_vars];
280        double initial_solution_lb = fp_solver_interface->getInfinity() * fp_solver_interface->getObjSense();
281
282        if (FPRandomInitialSolution == true)
283          {
284            DEBUGPR(10, ucout << "Starting search from random solution"
285                    << std::endl << Flush);
286
287            RNG * rng = pebbl::gRandomRNG();
288
289            for (int i = 0; i < num_vars; ++i)
290              {
291                if (lp->isContinuous(i) == true)
292                  {
293                    initial_solution[i] = 0.0; // isn't manipulated by FP anyway
294                  }
295                else
296                  {
297                    DUniform<int> sampler(rng,
298                                          int(rint(lp->getColLower(i))),
299                                          int(rint(lp->getColUpper(i))));
300                    initial_solution[i] = sampler();
301                  }
302              }
303          }
304        else
305          {
306            DEBUGPR(10, ucout << "Starting search from LP relaxation" << std::endl << Flush);
307
308            const double * x = fp_solver_interface->getColSolution();
309
310            DEBUGPR(100, prettyPrint("LP relaxation: ", x, num_vars, ucout));
311            DEBUGPR(100, ucout << "Objective=" << fp_solver_interface->getObjValue() << std::endl;);
312
313            initial_solution_lb = fp_solver_interface->getObjValue();
314           
315            // the first step in the feasibility pump is to solve the
316            // LP relaxation, in order to get the initial x^{*} vector.
317            // fortunately, we already have that as a side-effect of
318            // MILPNode. actually, any solution is OK - the idea behind
319            // using the LP relaxation is to have a prayer at finding
320            // a low-cost feasible solution.
321
322            for (int i = 0; i < num_vars; ++i)
323              {
324                // if a variable is continuous, leave it alone (it won't
325                // be examined anyway).  otherwise, round it to the
326                // nearest integer value.  NOTE: using lp instead of
327                // fp_solver_interface because there aren't any integers
328                // in the copy for some reason...
329                if (lp->isContinuous(i) == true)
330                  {
331                    initial_solution[i] = x[i]; // isn't used anyway
332                  }
333                else
334                  {
335                    initial_solution[i] = rint(x[i]);
336                  }
337              }
338          }
339       
340        // the feasibility pump objective is always to minimize distance.
341        fp_solver_interface->setObjSense(pebbl::minimize);
342
343        // add a row to the cloned solver interface to represent
344        // the objective cut. at the root node, this starts out at
345        // infinity, and is incrementally modified as better and
346        // better solutions are found - assuming incumbent improvement
347        // is enabled.
348        // NOTE: does milpNode or the base heuristic already add the objective cut?
349        const double * obj_coeffs = fp_solver_interface->getObjCoefficients();
350        CoinPackedVector objective_cut(num_vars, obj_coeffs);
351        fp_solver_interface->addRow(objective_cut, -(fp_solver_interface->getInfinity()), fp_solver_interface->getInfinity());
352        // the new row is added as the last index (num_rows - 1).
353        int num_rows(fp_solver_interface->getNumRows());
354
355        // call the core feasibility pump routine.
356        double * resulting_solution = new double[num_vars];
357        int iterations_this_trial(0);
358        double runtime_this_trial(0.0);
359        bool result = feasibilityPumpCore(fp_solver_interface,
360                                          initial_solution,
361                                          initial_solution_lb,
362                                          allowable_iterations
363                                          - cumulative_iterations,
364                                          FPRunTime
365                                          - (CPUSeconds() - start_time),
366                                          iterations_this_trial,
367                                          runtime_this_trial,
368                                          resulting_solution);
369
370        cumulative_iterations += iterations_this_trial;
371
372        if (result == true)
373          {
374            DEBUGPR(50, ucout << "Updating incumbent with feasibility pump"
375                    " solution" << std::endl << Flush);
376
377            // translate solution into utilib world.
378            DoubleVector heuristic_solution(num_vars);
379            for (int i = 0; i < num_vars; ++i)
380              {
381                heuristic_solution[i] = resulting_solution[i];
382              }
383
384            double this_objective = lp->evalVector(heuristic_solution);
385
386            if (((lp->getObjSense() == pebbl::minimize) &&
387                 (this_objective < best_objective)) ||
388                ((lp->getObjSense() == pebbl::maximize) &&
389                 (this_objective > best_objective)))
390              {
391                best_objective = this_objective;
392                best_solution = heuristic_solution;
393              }
394          }
395        else
396          {
397            break;
398          }
399
400        // clean up.
401        delete [] initial_solution;
402        delete [] resulting_solution;
403        delete fp_solver_interface;
404      }
405
406    MILPSolution * result = 0;
407
408    if (((lp->getObjSense() == pebbl::minimize) &&
409         (best_objective != std::numeric_limits<double>::max())) ||
410        ((lp->getObjSense() == pebbl::maximize) &&
411         (best_objective != std::numeric_limits<double>::min())))
412      {
413        DEBUGPR(10, ucout << "Objective value of best incumbent="
414                << best_objective << ", ");
415        DEBUGPR(10, ucout << "Meta-iterations=" << executed_iterations << ", ");
416        DEBUGPR(10, ucout << "Cumulative iterations=" << cumulative_iterations
417                << ", ");
418        DEBUGPR(10, ucout << "Cumulative CPU time=" << CPUSeconds() - start_time
419                << " seconds" << std::endl << Flush);
420
421        // NOTE: the heuristic solution is unlikely to be exactly integral
422        //       in the 'integer' variables. may cause problems when in
423        //       acro validation mode.
424
425        result = new MILPSolution(best_solution, mGlobal, lp->evalVector(best_solution));
426      }
427    else
428      {
429        DEBUGPR(10, ucout << "No feasible incumbent was found" << std::endl << Flush);
430      }
431
432    // revert the LP back to where it was ("solve" mode).
433    lp->setSolving();
434
435    DEBUGPR(100, ucout << "-----------------------------------------------"  << std::endl << Flush);
436
437
438    return result;
439  }
440
441  bool feasibilityPump::deterministicRounding(unsigned int & num_rounded_up,
442                                              unsigned int & num_rounded_down,
443                                              int num_vars,
444                                              const double * new_solution,
445                                              double * x_tilde)const
446  {
447    bool any_inverted = false;
448    num_rounded_down = 0;
449    num_rounded_up = 0;
450    for (int j = 0; j < num_vars; ++j)
451      {
452        if (lp->isContinuous(j) == false)
453          {
454            if (bGlobal->isInteger(new_solution[j]) == false)
455              {
456                double rounded_value = rint(new_solution[j]);
457                if (fabs(rounded_value - x_tilde[j]) > 1e-5) // enough of a threshold, given comparison of integers.
458                  {
459                    any_inverted = true;
460                    x_tilde[j] = rounded_value;
461                   
462                    // the following is simply for statistics.
463                    if (rounded_value <= new_solution[j])
464                      {
465                        ++num_rounded_down;
466                      }
467                    else
468                      {
469                        ++num_rounded_up;
470                      }
471
472                    // spit out the actual rounding information for debug/analysis purposes.
473                    if (verbosity(50) == true)
474                      {
475                        std::cout << "Rounded value=" << new_solution[j] << " to integral value=" << rounded_value << std::endl;
476                      }
477                  }
478              }
479          }
480      }
481
482    return any_inverted;
483  }
484
485  bool feasibilityPump::randomizedRounding(unsigned int & num_rounded_up,
486                                           unsigned int & num_rounded_down,
487                                           int num_vars,
488                                           const double * new_solution,
489                                           double * x_tilde)const
490  {
491    bool any_inverted = false;
492    num_rounded_down = 0;
493    num_rounded_up = 0;
494
495    RNG * rng = pebbl::gRandomRNG();
496    Uniform sampler(rng, 0.0, 1.0);
497
498    for (int j = 0; j < num_vars; ++j)
499      {
500        if (lp->isContinuous(j) == false)
501          {
502            if (bGlobal->isInteger(new_solution[j]) == false)
503              {
504                // if you're <= 0.33 or >= 0.67 in the fractional component,
505                // perform deterministic rounding. otherweise, round with
506                // equal probability in one or the other direction.
507                double fractional_component(new_solution[j] - floor(new_solution[j]));
508                double rounded_value;
509                if (fractional_component < 0.33)
510                  {
511                    rounded_value = floor(new_solution[j]);
512                  }
513                else if (fractional_component > 0.67)
514                  {
515                    rounded_value = ceil(new_solution[j]);                 
516                  }
517                else
518                  {
519                    if (sampler() <= 0.5)
520                      {
521                        rounded_value = floor(new_solution[j]);
522                      }
523                    else
524                      {
525                        rounded_value = ceil(new_solution[j]);
526                      }
527                  }
528
529                if (fabs(rounded_value - x_tilde[j]) > 1e-5) // enough of a threshold, given comparison of integers.
530                  {
531                    any_inverted = true;
532                    x_tilde[j] = rounded_value;
533                   
534                    // the following is simply for statistics.
535                    if (fabs(rounded_value - floor(new_solution[j])) < 1e-5)
536                      {
537                        ++num_rounded_down;
538                      }
539                    else
540                      {
541                        ++num_rounded_up;
542                      }
543
544                    // spit out the actual rounding information for debug/analysis purposes.
545                    if (verbosity(50) == true)
546                      {
547                        std::cout << "Rounded value=" << new_solution[j] << " to integral value=" << rounded_value << std::endl;
548                      }
549                  }
550              }
551          }
552      }
553
554    return any_inverted;
555  }
556
557  bool feasibilityPump::feasibilityPumpCore(OsiSolverInterface * fp_solver_interface,
558                                            const double * initial_solution,
559                                            double initial_solution_lb,
560                                            int allocated_iterations,
561                                            double allocated_runtime,
562                                            int & consumed_iterations,
563                                            double & consumed_runtime,
564                                            double * incumbent_solution)
565  {
566    assert(allocated_runtime >= 0.0);
567    assert(allocated_iterations >= 0);
568
569    double start_time = CPUSeconds();
570    int num_vars = fp_solver_interface->getNumCols();
571    int num_int_vars = problem->numBinaryVars();
572
573    consumed_iterations = 0;
574    consumed_runtime = 0.0;
575
576    if (FPRoundingType == 0)
577      {
578        DEBUGPR(10, ucout << "Deterministic rounding enabled for feasibility pump" << std::endl << Flush);
579      }
580    else // (FPRoundingType == 1)
581      {
582        DEBUGPR(10, ucout << "Randomized rounding enabled for feasibility pump" << std::endl << Flush);
583      }
584
585    // initialize the hash weights, which are used to detect cycles in the solution vector x tilde.
586    std::vector<int> hash_weights(num_vars);
587    for (int i = 0;i < num_vars; ++i)
588      {
589        RNG * rng = pebbl::gRandomRNG();
590        DUniform<int> sampler(rng, 1, std::numeric_limits<int>::max());
591        hash_weights[i] = sampler();
592      }
593   
594    const int CYCLE_HISTORY_LENGTH = 100; // could be promoted to a parameter at some point - should be linked to "R".
595    std::list<double> x_tilde_hash_history;
596
597    // the initial x~ is just the input solution.
598    double * x_tilde = new double[num_vars];
599    double hash_value(0.0);
600    for (int i = 0; i < num_vars; ++i)
601      {
602        x_tilde[i] = initial_solution[i];
603        hash_value += x_tilde[i] * hash_weights[i];
604      }
605
606    x_tilde_hash_history.push_back(hash_value);
607
608    // perform the specified number of feasibility pump iterations.
609    bool incumbent_found(false);
610    bool time_limited(false); // upon exit, is it due to time limiation?
611    bool iteration_limited(false); // upon exit, is it due to iteration limitation?
612    int iteration_feasible_found(std::numeric_limits<int>::min());
613
614    // the following is for reporting purposes only, and won't be triggered
615    // without some real weird inputs.
616    if (allocated_iterations == 0)
617      {
618        iteration_limited = true;
619      }
620
621    // in an attempt to accelerate solves...
622    CoinWarmStart * warm_start = 0;
623
624    double lp_relaxation_objective(fp_solver_interface->getInfinity() * fp_solver_interface->getObjSense());
625    double best_incumbent_objective(fp_solver_interface->getInfinity() * fp_solver_interface->getObjSense());
626
627    for (int current_iteration = 1; current_iteration <= allocated_iterations; ++current_iteration)
628      {
629        DEBUGPR(100, ucout << "Executing feasibility pump iteration=" << std::setw(5) << current_iteration << std::endl << Flush);
630
631        ++consumed_iterations;
632
633        DEBUGPR(100, prettyPrint("Reference x~ solution", x_tilde, num_vars, ucout));
634
635        // form the modified objective function for the LP, based on x~
636        unsigned base_offset(0);
637
638        for (int j = 0; j < num_vars; ++j)
639          {
640            // again, the lp instance has the integer information.
641            if (lp->isContinuous(j) == true)
642              {
643                fp_solver_interface->setObjCoeff(j, 0.0);
644              }
645            else // for now, we assume all integers are binary
646              {
647                if (x_tilde[j] == 0.0) // the rounded variables should be exactly equal to 0 (which is representable), so tolerances shouldn't be an issue.
648                  {
649                    fp_solver_interface->setObjCoeff(j, 1.0);
650                  }
651                else
652                  {
653                    base_offset++;
654                    fp_solver_interface->setObjCoeff(j, -1.0);
655                  }
656              }
657          }
658
659        const double * coeffs = fp_solver_interface->getObjCoefficients();
660        DEBUGPR(100, prettyPrint("New objective function coefficients in lp*", coeffs, num_vars, ucout));
661
662        // solve LP*
663        DEBUGPR(100, ucout << "Solving LP*" << std::endl << Flush);
664
665        // NOTE: I had originally tried an "initialSolve()" call on the first iteration,
666        //       and "resolve()" calls on subsequent iterations. this yielded some really
667        //       strange behaviors on some instances, including mas74 and opt1217. there
668        //       may be some confusion over what the term "problem modification" means
669        //       in OSI land, as the limited documentation indicates you can call "resolve"
670        //       after a problem modification.
671
672        if (warm_start != 0)
673          {
674            if (fp_solver_interface->setWarmStart(warm_start) == false)
675              {
676                DEBUGPR(10, ucout << "WARNING: COIN rejected the warm start object" << std::endl << Flush);
677              }
678          }
679
680        fp_solver_interface->messageHandler()->setLogLevel(0);
681        fp_solver_interface->initialSolve();
682
683        // save for the next round.
684        delete warm_start;
685        warm_start = fp_solver_interface->getWarmStart(); // ownership is trasnferred
686
687        DEBUGPR(100, ucout << "LP* objective value=" << fp_solver_interface->getObjValue() << std::endl << Flush);
688        DEBUGPR(100, ucout << "LP* base offset=" << base_offset << std::endl << Flush);
689        const double * new_solution = fp_solver_interface->getColSolution();
690        DEBUGPR(100, prettyPrint("LP* solution", new_solution, num_vars, ucout));
691
692        // compute the L1-norm distance between the new x* and the old x~
693        double dist(0.0);
694        for (int j = 0; j < num_vars; ++j)
695          {
696            if (lp->isContinuous(j) == false)
697              {
698                dist += fabs(x_tilde[j]-new_solution[j]);
699              }
700          }
701
702        DEBUGPR(10, ucout << "Iteration=" << std::setw(5) << current_iteration << ": Distance=" << std::setw(10) << dist << " " << ": Simplex iterations=" << fp_solver_interface->getIterationCount() <<  Flush);
703
704        // NOTE: The solution can still be integer-feasible, although the
705        //       distance is > 0. weird, but true. due to tolerances.
706
707        // see if the resulting solution is integer.
708        bool integer_feasible(true); // until proven otherwise
709        for (int j = 0; (integer_feasible == true) && (j < num_vars); ++j)
710          {
711            if (lp->isContinuous(j) == false)
712              {
713                // integer tolerances are found in the pebblBase
714                // class, as is this static utility method.
715                if (bGlobal->isInteger(new_solution[j]) == false)
716                  {
717                    integer_feasible = false;
718                  }
719              }
720          }
721
722        // if the solution is integer feasible, then update the best incumbent
723        // and update the objective cut accordingly (assuming continuation is enabled).
724        if (integer_feasible == true)
725          {
726            incumbent_found = true;
727            iteration_feasible_found = current_iteration;
728
729            // NOTE: this was the original - now keep going after modifying the objective function!
730            //      break;
731
732            const double * new_solution = fp_solver_interface->getColSolution();
733
734            DoubleVector heuristic_solution(num_vars);
735            for (int j = 0; j < num_vars; ++j)
736              {
737                heuristic_solution[j] = new_solution[j];
738                incumbent_solution[j] = new_solution[j]; // save the best-so-far.
739              }
740
741            double solution_objective = lp->evalVector(heuristic_solution);
742            best_incumbent_objective = solution_objective;
743
744            DEBUGPR(10, ucout << ", >>>>solution is integer feasible - objective=" << solution_objective << std::endl << Flush);
745
746            x_tilde_hash_history.clear();
747
748            // only continue search if the improve-incumbent flag is enabled -
749            // otherwise, bail once you have a feasible incumbent.
750            if (FPImproveIncumbent == false)
751              {
752                DEBUGPR(10, ucout << "Incumbent improvement is disabled - terminating search" << std::endl << Flush);
753                break;
754              }
755
756            // TBD - still need to mess with the objective sense - it's currently hardwired for minimization.
757            double absolute_gap = solution_objective - initial_solution_lb;
758            double new_objective = solution_objective - FPImprovementAlpha * absolute_gap;
759            int num_rows = fp_solver_interface->getNumRows();
760            fp_solver_interface->setRowUpper(num_rows-1, new_objective);
761
762            DEBUGPR(10, ucout << "Objective upper bound reset to " << new_objective << std::endl << Flush);
763          }
764        else
765          {
766            // see if there is a cycle in the x tilde values.
767            bool cycle_detected(false);
768            unsigned int cycle_length(0);
769            std::list<double>::reverse_iterator history_iterator(x_tilde_hash_history.rbegin());
770            double current_hash_value(*history_iterator);
771            for (++history_iterator; history_iterator != x_tilde_hash_history.rend(); ++history_iterator)
772              {
773                const int MINIMUM_CYCLE_LENGTH = 5; // random moves can get you out of simple cases, so let it go a bit.
774                ++cycle_length;
775                if ((cycle_length >= MINIMUM_CYCLE_LENGTH) &&
776                    (fabs((*history_iterator) - current_hash_value) < 1e-5))
777                  {
778                    DEBUGPR(10, ucout << ", Cycle detected in x tilde vector, cycle length=" << cycle_length << Flush);
779                    cycle_detected = true;
780                    break;
781                  }
782              }
783
784            // NOTE: there is a hack described in the original feasibility pump
785            //       paper, in which a random move is applied every FPPerturbationInterval
786            // iterations to generate x~, instead of using x*.
787
788            if (((current_iteration % FPPerturbationInterval) == 0) ||
789                (cycle_detected == true))
790              {
791                DEBUGPR(10, ucout << ", Applying perturbation" << std::endl << Flush);
792     
793                RNG * rng = pebbl::gRandomRNG();
794                Uniform sampler(rng, -0.3, 0.7);
795     
796                for (int j = 0; j < num_vars; ++j)
797                  {
798                    if (lp->isContinuous(j) == false)
799                      {
800                        double sample = sampler();
801                        double result = fabs(new_solution[j] - x_tilde[j]) + std::max(sample,0.0);
802                 
803                        if (result >= 0.5)
804                          {
805                            x_tilde[j] = (1.0 - x_tilde[j]);
806                          }
807                      }
808                  }
809
810                x_tilde_hash_history.clear();
811              }   
812            else
813              {
814                // construct the next x tilde via any number of mechanisms, to see if rounding
815                // of any x* variables yields new values for the corresponding x~ variables.
816
817                bool any_inverted(false);
818                unsigned int num_rounded_up(0);
819                unsigned int num_rounded_down(0);
820
821                if (FPRoundingType == 0)
822                  {
823                    any_inverted = deterministicRounding(num_rounded_up,
824                                                         num_rounded_down,
825                                                         num_vars,
826                                                         new_solution,
827                                                         x_tilde);
828                  }
829                else // (FPRoundingType == 1)
830                  {
831                    any_inverted = randomizedRounding(num_rounded_up,
832                                                      num_rounded_down,
833                                                      num_vars,
834                                                      new_solution,
835                                                      x_tilde);
836                  }
837
838                if (any_inverted == true)
839                  {
840                    DEBUGPR(10, ucout << ", Variables were inverted - " << num_rounded_down << " rounded down, " << num_rounded_up << " rounded up" << std::endl << Flush);
841                  }
842                else
843                  {
844                    // if the new x* rounding failed to change x~,
845                    // flip some number of variables with the highest
846                    // value of abs(x* - x~).
847                    DEBUGPR(100, ucout << "Mean flip count=" << FPMeanFlipCount << std::endl << Flush);
848
849                    // compute |x* - x~| for each variable and sort in descending order.
850                    std::vector<std::pair<double,int> > distances;
851                    distances.reserve(num_vars);
852                    for (int j = 0; j < num_vars; ++j)
853                      {
854                        if (lp->isContinuous(j) == false)
855                          {
856                            double this_distance = fabs(new_solution[j]-x_tilde[j]);
857                            distances.push_back(std::make_pair(this_distance,j));
858                          }
859                      }
860
861                    std::sort(distances.begin(),distances.end(),compareDists);
862
863                    //              std::cout << "Distance/flip info: ";
864                    //              for (size_t i = 0; i < distances.size(); ++i)
865                    //                {
866                    //                  std::cout << "(" << distances[i].first << "," << distances[i].second << ") ";
867                    //                }
868                    //              std::cout << std::endl;
869
870                    // flip the rand(fpFlipCount/2,3*fpFlipCount/2) variables with the largest distance
871                    RNG * rng = pebbl::gRandomRNG();
872                    DUniform<int> sampler(rng,
873                                          std::max(int(floor(double(FPMeanFlipCount)/2.0)),1),
874                                          int(ceil(3.0*double(FPMeanFlipCount)/2.0)));
875
876                    int num_to_flip = sampler();
877                    // can't flip more variables than there are!
878                    if (num_to_flip > num_int_vars)
879                      {
880                        num_to_flip = num_int_vars;
881                      }
882                   
883                    DEBUGPR(100, ucout << "Actual number of flips to perform=" << num_to_flip << std::endl << Flush);
884
885                    // normalize flip count to the maximal number of binary variables.
886                    // NOTE: a lot of these are probably the same, and we should randomize within the blocks.
887                    num_to_flip = std::min(num_to_flip,int(distances.size()));
888
889                    for (int j = 0; j < num_to_flip; ++j)
890                      {
891                        int this_var = distances[j].second;
892                        //                      std::cout << "This distance=" << distances[j].first << std::endl;
893                        x_tilde[this_var] = 1.0 - x_tilde[this_var];
894                      }
895
896                    DEBUGPR(10, ucout << ", Randomly flipped " << num_to_flip << " variables" << std::endl << Flush);
897
898                    DEBUGPR(100, prettyPrint("X~ solution after randomized flipping", x_tilde, num_vars, ucout));
899                  }
900              }
901
902            // update x tilde history for cycle detection.
903            hash_value = 0.0;
904            for (int i = 0; i < num_vars; ++i)
905              {
906                hash_value += x_tilde[i] * hash_weights[i];
907              }
908
909            x_tilde_hash_history.push_back(hash_value);
910           
911            while (int(x_tilde_hash_history.size()) > CYCLE_HISTORY_LENGTH)
912              {
913                x_tilde_hash_history.pop_front();
914              }
915          }
916
917        double elapsed_time = CPUSeconds() - start_time;
918        if (elapsed_time >= allocated_runtime)
919          {
920            time_limited = true;
921
922            break;
923          }
924
925        if (current_iteration == allocated_iterations)
926          {
927            DEBUGPR(10, ucout << "Feasibility pump terminating after exceeeding allocated number of iterations" << std::endl << Flush);
928
929            iteration_limited = true;
930            // the loop will terminate naturally at this point - no break is necessary.
931          }
932
933        DEBUGPR(50, ucout << std::endl << Flush);
934      }
935
936    double end_time = CPUSeconds();
937
938    if (incumbent_found == true)
939      {
940        DoubleVector solution_to_evaluate(num_vars);
941        for (int j = 0; j < num_vars; ++j)
942          {
943            solution_to_evaluate[j] = incumbent_solution[j];
944          }
945
946        double solution_objective = lp->evalVector(solution_to_evaluate);
947
948        DEBUGPR(10, ucout << "Feasibility pump found feasible integer solution: Iterations=" << iteration_feasible_found << ", ");
949        DEBUGPR(10, ucout << "Objective value=" << solution_objective << ", ");
950        DEBUGPR(10, ucout << "CPU time=" << end_time - start_time << " seconds" << std::endl << Flush);
951      }
952    else
953      {
954        DEBUGPR(10, ucout << std::endl << "Feasibility pump failed to find feasible integer solution: ");
955        if (iteration_limited == true)
956          {
957            DEBUGPR(10, ucout << "Iteration limit=" << allocated_iterations << " exceeded");
958          }
959        else // time_limited == true
960          {
961            DEBUGPR(10, ucout << "Time limit=" << allocated_runtime << " seconds exceeded");
962          }
963        DEBUGPR(10, ucout << std::endl << Flush);
964      }
965
966    delete warm_start;
967
968    delete [] x_tilde;
969
970    consumed_runtime = CPUSeconds() - start_time;
971
972    return (incumbent_found == true);
973  }
974
975 
976} // namespace pico
Note: See TracBrowser for help on using the repository browser.