itm-common.hpp 57.1 KB
Newer Older
Gauthier Quesnel's avatar
2019    
Gauthier Quesnel committed
1
/* Copyright (C) 2016-2019 INRA
Gauthier Quesnel's avatar
Gauthier Quesnel committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

#ifndef ORG_VLEPROJECT_BARYONYX_SOLVER_ITM_COMMON_HPP
#define ORG_VLEPROJECT_BARYONYX_SOLVER_ITM_COMMON_HPP

26
27
#include <baryonyx/core>

28
#include "bit-array.hpp"
29
#include "debug.hpp"
30
#include "private.hpp"
31
#include "problem.hpp"
32
#include "sparse-matrix.hpp"
33
#include "utils.hpp"
34
35
36

#include <algorithm>
#include <chrono>
37
#include <fstream>
38
#include <iterator>
Gauthier Quesnel's avatar
Gauthier Quesnel committed
39
#include <numeric>
40
#include <random>
41
#include <thread>
42
#include <tuple>
43
#include <vector>
44

45
46
#include <fmt/ostream.h>

Gauthier Quesnel's avatar
Gauthier Quesnel committed
47
48
49
namespace baryonyx {
namespace itm {

50
51
using random_engine = std::default_random_engine;

52
53
54
55
56
57
58
struct maximize_tag
{};

struct minimize_tag
{};

struct merged_constraint
59
{
60
61
62
63
64
65
66
67
    merged_constraint(std::vector<function_element> elements_,
                      int min_,
                      int max_,
                      int id_)
      : elements(std::move(elements_))
      , min(min_)
      , max(max_)
      , id(id_)
68
    {}
69
70
71
72
73

    std::vector<function_element> elements;
    int min;
    int max;
    int id;
74
75
};

76
77
std::vector<merged_constraint>
make_merged_constraints(const context_ptr& ctx, const problem& pb);
78

79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
template<typename Solver, typename Xtype>
bool
is_valid_constraint(const Solver& slv, int k, const Xtype& x)
{
    typename sparse_matrix<int>::const_row_iterator it, et;

    std::tie(it, et) = slv.ap.row(k);
    int v = 0;

    for (; it != et; ++it)
        v += slv.factor(it->value) * x[it->column];

    return slv.bound_min(k) <= v && v <= slv.bound_max(k);
}

template<typename Solver, typename Xtype>
bool
is_valid_solution(const Solver& slv, const Xtype& x)
{
    for (int k = 0; k != slv.m; ++k)
        if (!is_valid_constraint(slv, k, x))
            return false;

    return true;
}

template<typename Solver, typename Xtype>
int
compute_violated_constraints(const Solver& slv,
                             const Xtype& x,
109
                             std::vector<int>& out)
110
111
112
113
114
115
116
117
118
119
{
    out.clear();

    for (int k = 0; k != slv.m; ++k)
        if (!is_valid_constraint(slv, k, x))
            out.emplace_back(k);

    return length(out);
}

120
template<typename iteratorT>
121
inline void
122
123
124
random_shuffle_unique(iteratorT begin,
                      iteratorT end,
                      random_engine& rng) noexcept
125
126
127
128
129
130
131
132
133
134
135
136
{
    auto ret = begin++;
    for (; begin != end; ++begin) {
        if (ret->value != begin->value) {
            std::shuffle(ret, begin, rng);
            ret = begin;
        }
    }

    std::shuffle(ret, begin, rng);
}

137
template<typename Mode, typename Iterator>
138
inline void
139
calculator_sort(Iterator begin, Iterator end, random_engine& rng)
140
141
142
{
    if (std::distance(begin, end) > 1) {
        std::sort(begin, end, [](const auto& lhs, const auto& rhs) {
143
            if constexpr (std::is_same_v<Mode, minimize_tag>)
144
145
146
                return lhs.value < rhs.value;
            else
                return rhs.value < lhs.value;
147
148
149
150
151
152
        });

        random_shuffle_unique(begin, end, rng);
    }
}

153
template<typename Mode, typename Float>
154
inline bool
155
stop_iterating(Float value, random_engine& rng) noexcept
156
157
158
159
160
161
{
    if (value == 0) {
        std::bernoulli_distribution d(0.5);
        return d(rng);
    }

162
163
164
165
    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return value > 0;
    else
        return value < 0;
166
167
}

168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
template<typename Mode, typename Float>
constexpr Float
bad_value() noexcept
{
    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return std::numeric_limits<Float>::infinity();
    else
        return -std::numeric_limits<Float>::infinity();
}

template<typename Mode, typename Float>
constexpr bool
is_bad_value(Float value) noexcept
{
    return value == bad_value<Mode, Float>();
}

185
template<typename Mode, typename Float>
186
inline bool
187
stop_iterating(Float value) noexcept
188
{
189
190
191
192
    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return value > 0;
    else
        return value < 0;
193
194
}

195
template<typename Mode, typename Float>
196
inline bool
197
is_better_solution(Float lhs, Float rhs) noexcept
198
{
199
200
201
202
    if constexpr (std::is_same_v<Mode, minimize_tag>)
        return lhs < rhs;
    else
        return lhs > rhs;
203
204
}

205
template<typename Mode, typename Float>
206
inline bool
207
init_x(Float cost, int value_if_cost_0) noexcept
208
{
209
210
211
    if constexpr (std::is_same_v<Mode, minimize_tag>) {
        if (cost < 0)
            return true;
212

213
214
215
216
217
        if (cost == 0)
            return value_if_cost_0;
    } else {
        if (cost > 0)
            return true;
218

219
220
221
        if (cost == 0)
            return value_if_cost_0;
    }
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249

    return false;
}

template<typename TimePoint>
double
compute_duration(const TimePoint& first, const TimePoint& last) noexcept
{
    namespace sc = std::chrono;

    auto diff = last - first;
    auto dc = sc::duration_cast<sc::duration<double, std::ratio<1>>>(diff);

    return dc.count();
}

inline std::size_t
compute_reduced_costs_vector_size(
  const std::vector<merged_constraint>& csts) noexcept
{
    std::size_t rsizemax = csts[0].elements.size();

    for (std::size_t i = 1, e = csts.size(); i != e; ++i)
        rsizemax = std::max(rsizemax, csts[i].elements.size());

    return rsizemax;
}

250
251
252
253
254
255
256
257
258
259
260
261
262
263
/// Test if the bit sign between two reals are differents.
template<typename Float>
inline bool
is_signbit_change(Float lhs, Float rhs) noexcept
{
    return std::signbit(lhs) != std::signbit(rhs);
}

/// The bkmin and bkmax constraint bounds are not equal and can be assigned to
/// -infinity or +infinity. We have to scan the r vector and search a value j
/// such as b(0, k) <= Sum A(k, R[j]) < b(1, k).
///
/// @return True if the sign of the pi vector change during the affect
///     operation.
264
template<typename Solver, typename Xtype, typename Iterator, typename Float>
265
bool
266
267
268
269
270
271
272
273
274
275
276
277
278
affect(Solver& slv,
       Xtype& x,
       Iterator it,
       int k,
       int selected,
       int r_size,
       const Float kappa,
       const Float delta)
{
    constexpr Float one{ 1 };
    constexpr Float two{ 2 };
    constexpr Float middle{ (two + one) / two };

279
280
    const auto old_pi = slv.pi[k];

281
282
283
    auto d = delta;

    if (selected < 0) {
284
        // slv.pi[k] += slv.R[0].value / two;
285
286
287
288
289
        d += (kappa / (one - kappa)) * (slv.R[0].value / two);

        for (int i = 0; i != r_size; ++i) {
            auto var = it + slv.R[i].id;

290
            if (slv.R[i].is_negative_factor()) {
291
                x.set(var->column);
292
293
                slv.P[var->value] += d;
            } else {
294
                x.unset(var->column);
295
296
297
298
                slv.P[var->value] -= d;
            }
        }
    } else if (selected + 1 >= r_size) {
299
        // slv.pi[k] += slv.R[selected].value * middle;
300
301
302
303
304
        d += (kappa / (one - kappa)) * (slv.R[selected].value * middle);

        for (int i = 0; i != r_size; ++i) {
            auto var = it + slv.R[i].id;

305
            if (slv.R[i].is_negative_factor()) {
306
                x.unset(var->column);
307
308
                slv.P[var->value] -= d;
            } else {
309
                x.set(var->column);
310
311
312
313
314
315
316
317
318
319
320
321
322
                slv.P[var->value] += d;
            }
        }
    } else {
        slv.pi[k] +=
          ((slv.R[selected].value + slv.R[selected + 1].value) / two);
        d += (kappa / (one - kappa)) *
             (slv.R[selected + 1].value - slv.R[selected].value);

        int i = 0;
        for (; i <= selected; ++i) {
            auto var = it + slv.R[i].id;

323
            if (slv.R[i].is_negative_factor()) {
324
                x.unset(var->column);
325
326
                slv.P[var->value] -= d;
            } else {
327
                x.set(var->column);
328
329
330
331
332
333
334
                slv.P[var->value] += d;
            }
        }

        for (; i != r_size; ++i) {
            auto var = it + slv.R[i].id;

335
            if (slv.R[i].is_negative_factor()) {
336
                x.set(var->column);
337
338
                slv.P[var->value] += d;
            } else {
339
                x.unset(var->column);
340
341
342
343
344
345
346
                slv.P[var->value] -= d;
            }
        }
    }

    // TODO job: develops is_valid_constraint for all the solvers
    bx_expects(is_valid_constraint(slv, k, x));
347
348

    return is_signbit_change(old_pi, slv.pi[k]);
349
350
}

351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
namespace detail {

inline int
constraint(std::vector<int>::const_iterator it)
{
    return *it;
}

inline int
constraint(std::vector<int>::const_reverse_iterator it)
{
    return *it;
}

inline int
constraint(std::vector<std::pair<int, int>>::const_iterator it)
{
    return it->first;
}
}

template<typename Iterator>
int
constraint(Iterator it)
{
    return detail::constraint(it);
}

379
380
381
template<typename Float, typename Mode>
struct best_solution_recorder;

382
template<typename Solver, typename Float, typename Mode>
383
class solver_initializer
384
{
385
    std::bernoulli_distribution dist;
386
    std::bernoulli_distribution toss_up;
387
    bool use_cycle{ true };
388

389
    enum class single_automaton : std::uint8_t
390
    {
391
        init,
392
393
394
        improve_x_1,
        improve_x_2,
        improve_x_3
395
    };
396

397
    enum class cycle_automaton : std::uint8_t
398
    {
399
        bastert,
400
        bastert_crossover,
401
        pessimistic_solve,
402
        pessimistic_crossover,
403
        optimistic_solve,
404
        optimistic_crossover
405
406
    };

407
408
    std::array<unsigned, 4> cycle_counter = { 0, 0, 0, 0 };

409
410
411
412
413
414
415
    enum class init_automaton : std::uint8_t
    {
        random,
        current_best,
        best_recorder
    };

416
    static inline const std::string_view single_automaton_string[] =
417
      { "init", "improve_x_1", "improve_x_2", "improve_x_3" };
418
419

    static inline const std::string_view cycle_automaton_string[] = {
420
421
422
        "bastert",           "bastert_crossover",
        "pessimistic_solve", "pessimistic_solve_crossver",
        "optimistic_solve",  "optimistic_solve_crossover",
423
    };
424

425
426
    static inline const std::string_view init_automaton_string[] =
      { "random", "current_best", "best_recorder" };
Gauthier Quesnel's avatar
Gauthier Quesnel committed
427

428
429
    single_automaton single{ single_automaton::init };
    cycle_automaton cycle{ cycle_automaton::pessimistic_solve };
430
431
432
433
434
435
436
437
438
439
    init_automaton base{ init_automaton::random };

    void init_crossover(
      Solver& slv,
      bit_array& x,
      const best_solution_recorder<Float, Mode>& best_recorder)
    {
        to_log(stdout, 6u, "initializer: crossover {}\n", cycle_counter[3]);
        ++cycle_counter[3];

440
441
442
443
444
        const int min = 1;
        const int max = slv.n >= 2 ? slv.n - 1 : 1;

        std::uniform_int_distribution<int> uniform_dist(min, max);
        const int split_at = uniform_dist(slv.rng);
445
446
447
448

        // TODO Need to improve this part to avoid bit-per-bit copy using the
        // std::uint64_t underlying type.

449
450
451
452
453
454
455
456
        if (toss_up(slv.rng)) {
            if (toss_up(slv.rng)) {
                best_recorder.copy_any_solution(x, slv.rng, 0, split_at);
                best_recorder.copy_best_solution(x, split_at, slv.n);
            } else {
                best_recorder.copy_best_solution(x, 0, split_at);
                best_recorder.copy_any_solution(x, slv.rng, split_at, slv.n);
            }
457
        } else {
458
459
            best_recorder.copy_any_solution(x, slv.rng, 0, split_at);
            best_recorder.copy_any_solution(x, slv.rng, split_at, slv.n);
460
461
        }
    }
462

463
    void init_bastert(Solver& slv, bit_array& x) noexcept
464
    {
465
466
467
        to_log(stdout, 6u, "- initializer: bastert {}\n", cycle_counter[0]);
        ++cycle_counter[0];

468
469
470
        const int value_if_cost_0 = 1;

        for (int i = 0; i != slv.n; ++i) {
471
472
473
474
475
476
            if (dist(slv.rng)) {
                if (init_x<Mode>(slv.c(i, x), value_if_cost_0))
                    x.set(i);
                else
                    x.unset(i);
            }
477
        }
478
479
    }

480
    void init_pessimistic_solve(Solver& slv, bit_array& x) noexcept
481
    {
482
483
484
485
        to_log(
          stdout, 6u, "- initializer: pessimistic {}\n", cycle_counter[1]);
        ++cycle_counter[1];

486
487
488
        for (int k = 0; k != slv.m; ++k) {
            if (!dist(slv.rng))
                continue;
489

490
491
492
493
            auto [begin, end] = slv.ap.row(k);
            int r_size = 0;

            for (auto it = begin; it != end; ++it, ++r_size) {
Gauthier Quesnel's avatar
Gauthier Quesnel committed
494
                slv.R[r_size].value = slv.c(it->column, x);
495
496
497
                slv.R[r_size].id = r_size;
            }

498
499
            std::shuffle(slv.R.get(), slv.R.get() + r_size, slv.rng);

500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
            std::sort(slv.R.get(),
                      slv.R.get() + r_size,
                      [](const auto& lhs, const auto& rhs) {
                          if constexpr (std::is_same_v<Mode, minimize_tag>)
                              return lhs.value < rhs.value;
                          else
                              return rhs.value < lhs.value;
                      });

            int sum = 0;
            int best = -2;

            for (int i = -1; i < r_size; ++i) {
                if (slv.bound_min(k) <= sum && sum <= slv.bound_max(k)) {
                    best = i;
                    break;
                }

518
519
                auto var = begin + slv.R[i].id;
                sum += slv.factor(var->value);
520
521
522
523
524
            }

            int i = 0;
            for (; i <= best; ++i) {
                auto var = begin + slv.R[i].id;
525
                x.set(var->column);
526
527
528
529
            }

            for (; i != r_size; ++i) {
                auto var = begin + slv.R[i].id;
530
                x.unset(var->column);
531
532
533
534
            }
        }
    }

535
    void init_optimistic_solve(Solver& slv, bit_array& x) noexcept
536
    {
537
538
539
        to_log(stdout, 6u, "- initializer: optimistic {}\n", cycle_counter[2]);
        ++cycle_counter[2];

540
541
542
543
544
545
546
547
        for (int k = 0; k != slv.m; ++k) {
            if (!dist(slv.rng))
                continue;

            auto [begin, end] = slv.ap.row(k);
            int r_size = 0;

            for (auto it = begin; it != end; ++it, ++r_size) {
Gauthier Quesnel's avatar
Gauthier Quesnel committed
548
                slv.R[r_size].value = slv.c(it->column, x);
549
550
551
                slv.R[r_size].id = r_size;
            }

552
553
            std::shuffle(slv.R.get(), slv.R.get() + r_size, slv.rng);

554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
            std::sort(slv.R.get(),
                      slv.R.get() + r_size,
                      [](const auto& lhs, const auto& rhs) {
                          if constexpr (std::is_same_v<Mode, minimize_tag>)
                              return lhs.value < rhs.value;
                          else
                              return rhs.value < lhs.value;
                      });

            int sum = 0;
            int best = -2;

            for (int i = -1; i < r_size; ++i) {
                if (slv.bound_min(k) <= sum && sum <= slv.bound_max(k))
                    best = i;

                if (best != -2 && i + 1 < r_size &&
                    stop_iterating<Mode>(slv.R[i + 1].value))
                    break;

574
575
                auto var = begin + slv.R[i].id;
                sum += slv.factor(var->value);
576
577
578
579
580
            }

            int i = 0;
            for (; i <= best; ++i) {
                auto var = begin + slv.R[i].id;
581
                x.set(var->column);
582
583
584
585
            }

            for (; i != r_size; ++i) {
                auto var = begin + slv.R[i].id;
586
                x.unset(var->column);
587
588
589
590
            }
        }
    }

591
592
593
public:
    solver_initializer(Solver& slv,
                       solver_parameters::init_policy_type policy,
594
                       double init_policy_random,
595
                       double init_random)
596
597
      : dist(init_policy_random)
      , toss_up(init_random)
598
    {
599
        // x.clear();
600
        slv.reset();
601

602
        switch (policy) {
603
        case solver_parameters::init_policy_type::bastert:
604
            cycle = cycle_automaton::bastert;
605
606
607
608
609
            use_cycle = false;
            break;
        case solver_parameters::init_policy_type::pessimistic_solve:
            cycle = cycle_automaton::pessimistic_solve;
            use_cycle = false;
610
            break;
611
612
613
        case solver_parameters::init_policy_type::optimistic_solve:
            cycle = cycle_automaton::optimistic_solve;
            use_cycle = false;
614
            break;
615
616
617
        case solver_parameters::init_policy_type::cycle:
            cycle = cycle_automaton::pessimistic_solve;
            use_cycle = true;
618
            break;
Gauthier Quesnel's avatar
Gauthier Quesnel committed
619
        }
620
    }
Gauthier Quesnel's avatar
Gauthier Quesnel committed
621

622
    ~solver_initializer() noexcept
623
    {
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
        to_log(stdout,
               6u,
               "- solver-initializer status: bastert {} /"
               "pessimistic_solve {} / optimistic_solve {} / crossover {}\n",
               cycle_counter[0],
               cycle_counter[1],
               cycle_counter[2],
               cycle_counter[3]);
    }

    void next_state(bool is_a_solution) noexcept
    {
        if (is_a_solution) {
            single = single_automaton::init;
        } else {
            switch (single) {
            case single_automaton::init:
641
                single = single_automaton::improve_x_1;
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
                break;
            case single_automaton::improve_x_1:
                single = single_automaton::improve_x_2;
                break;
            case single_automaton::improve_x_2:
                single = single_automaton::improve_x_3;
                break;
            case single_automaton::improve_x_3:
                single = single_automaton::init;
            }
        }

        if (single == single_automaton::init) {
            if (use_cycle) {
                switch (cycle) {
                case cycle_automaton::bastert:
                    cycle = cycle_automaton::bastert_crossover;
                    break;
                case cycle_automaton::bastert_crossover:
                    cycle = cycle_automaton::pessimistic_solve;
                    break;
                case cycle_automaton::pessimistic_solve:
                    cycle = cycle_automaton::pessimistic_crossover;
                    break;
                case cycle_automaton::pessimistic_crossover:
                    cycle = cycle_automaton::optimistic_solve;
                    break;
                case cycle_automaton::optimistic_solve:
                    cycle = cycle_automaton::optimistic_crossover;
                    break;
                case cycle_automaton::optimistic_crossover:
                    cycle = cycle_automaton::bastert;
                    break;
                }
676
            } else {
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
                switch (cycle) {
                case cycle_automaton::bastert:
                    cycle = cycle_automaton::bastert_crossover;
                    break;
                case cycle_automaton::bastert_crossover:
                    cycle = cycle_automaton::bastert;
                    break;
                case cycle_automaton::pessimistic_solve:
                    cycle = cycle_automaton::pessimistic_crossover;
                    break;
                case cycle_automaton::pessimistic_crossover:
                    cycle = cycle_automaton::pessimistic_solve;
                    break;
                case cycle_automaton::optimistic_solve:
                    cycle = cycle_automaton::optimistic_crossover;
                    break;
                case cycle_automaton::optimistic_crossover:
                    cycle = cycle_automaton::optimistic_solve;
                    break;
696
697
                }
            }
698
        }
699
700
    }

701
    void init(Solver& slv, bit_array& x)
702
    {
703
        // x.clear();
704
        slv.reset();
705

706
        to_log(stdout,
707
708
               6u,
               "- initializer: init method {} in mode {} - init {}\n",
709
               cycle_automaton_string[static_cast<int>(cycle)],
710
               init_automaton_string[static_cast<int>(base)],
711
               single_automaton_string[static_cast<int>(single)]);
712

713
714
715
716
717
718
719
720
721
        for (int i = 0; i != slv.n; ++i)
            if (toss_up(slv.rng))
                x.set(i);
            else
                x.unset(i);

        switch (cycle) {
        case cycle_automaton::bastert:
            init_bastert(slv, x);
722
            break;
723
724
725
726
727
728
729
730
731
732
        case cycle_automaton::pessimistic_solve:
            init_pessimistic_solve(slv, x);
            break;
        case cycle_automaton::optimistic_solve:
            init_optimistic_solve(slv, x);
            break;
        case cycle_automaton::bastert_crossover:
        case cycle_automaton::pessimistic_crossover:
        case cycle_automaton::optimistic_crossover:
            bx_reach();
733
            break;
734
735
736
        }
    }

737
738
739
740
741
    Float reinit(Solver& slv,
                 bit_array& x,
                 bool is_a_solution,
                 const best_solution_recorder<Float, Mode>& best_recorder,
                 double kappa_min,
742
                 double kappa_max)
743
744
745
746
747
    {
        // x.clear();
        slv.reset();

        to_log(stdout,
748
749
750
               6u,
               "- initializer: re-init method {} in mode {} - init {} - "
               "is-a-solution {}\n",
751
               cycle_automaton_string[static_cast<int>(cycle)],
752
753
754
755
756
               init_automaton_string[static_cast<int>(base)],
               single_automaton_string[static_cast<int>(single)],
               is_a_solution);

        next_state(is_a_solution);
757

758
759
        switch (single) {
        case single_automaton::init:
760
            switch (base) {
761
762
763
764
765
766
767
            case init_automaton::random:
                for (int i = 0; i != slv.n; ++i)
                    if (toss_up(slv.rng))
                        x.set(i);
                    else
                        x.unset(i);

768
                base = init_automaton::current_best;
769
770
                break;
            case init_automaton::current_best:
771
772
773
774
775
776
777
778
                if (!best_recorder.copy_any_solution(x, slv.rng)) {
                    for (int i = 0; i != slv.n; ++i)
                        if (toss_up(slv.rng))
                            x.set(i);
                        else
                            x.unset(i);
                }

779
                base = init_automaton::best_recorder;
780
781
                break;
            case init_automaton::best_recorder:
782
783
784
785
786
787
788
789
                if (!best_recorder.copy_best_solution(x)) {
                    for (int i = 0; i != slv.n; ++i)
                        if (toss_up(slv.rng))
                            x.set(i);
                        else
                            x.unset(i);
                }

790
                base = init_automaton::random;
791
792
793
794
795
796
797
798
799
800
801
802
803
                break;
            }

            switch (cycle) {
            case cycle_automaton::bastert:
                init_bastert(slv, x);
                break;
            case cycle_automaton::pessimistic_solve:
                init_pessimistic_solve(slv, x);
                break;
            case cycle_automaton::optimistic_solve:
                init_optimistic_solve(slv, x);
                break;
804
805
806
            case cycle_automaton::bastert_crossover:
            case cycle_automaton::pessimistic_crossover:
            case cycle_automaton::optimistic_crossover:
807
                init_crossover(slv, x, best_recorder);
808
809
                break;
            }
810
811
            break;
        case single_automaton::improve_x_1:
812
            kappa_min += (kappa_max - kappa_min) / 20;
813
814
            break;
        case single_automaton::improve_x_2:
815
            kappa_min += (kappa_max - kappa_min) / 10;
816
817
            break;
        case single_automaton::improve_x_3:
818
            kappa_min += (kappa_max - kappa_min) / 5;
819
            break;
820
        }
821
822

        return static_cast<Float>(kappa_min);
823
824
    }
};
825

826
/**
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
 * Compute a problem lower or upper bounds based on Lagrangian multipliers
 * (valid if there are equality constraints only?)
 */
template<typename floatingpointT, typename modeT>
struct bounds_printer
{
    floatingpointT bestlb;
    floatingpointT bestub;
    floatingpointT max_cost;

    bounds_printer(floatingpointT max_cost_init)
      : bestlb(std::numeric_limits<floatingpointT>::lowest())
      , bestub(std::numeric_limits<floatingpointT>::max())
      , max_cost(max_cost_init)
    {}

    template<typename SolverT>
    floatingpointT init_bound(const SolverT& slv)
    {
        floatingpointT b{ 0 };

        for (auto c = 0; c != slv.m; ++c)
            b += slv.pi[c] * static_cast<floatingpointT>(slv.bound_init(c));

        return b;
    }

    template<typename SolverT>
    floatingpointT add_bound(const SolverT& slv,
                             int j,
                             floatingpointT sum_a_pi,
                             minimize_tag)
    {
        if (slv.c[j] - sum_a_pi < 0)
            return slv.c[j] - sum_a_pi;

        return { 0 };
    }

    template<typename SolverT>
    floatingpointT add_bound(const SolverT& slv,
                             int j,
                             floatingpointT sum_a_pi,
                             maximize_tag)
    {
        if (slv.c[j] - sum_a_pi > 0)
            return slv.c[j] - sum_a_pi;

        return { 0 };
    }

    void print_bound(const context_ptr& ctx,
                     floatingpointT lower_bound,
                     floatingpointT upper_bound,
                     minimize_tag)
    {
        bool better_gap = (lower_bound > bestlb || upper_bound < bestub);

        if (upper_bound < bestub)
            bestub = upper_bound;

        if (lower_bound > bestlb)
            bestlb = lower_bound;

        if (better_gap) {
            if (bestub == static_cast<floatingpointT>(0.0))
                info(ctx, "  - Lower bound: {}   (gap: 0%)\n", bestlb);
            else
                info(ctx,
                     "  - Lower bound: {}   (gap: {}%)\n",
                     bestlb,
                     static_cast<floatingpointT>(100.) * (bestub - bestlb) /
                       bestub);
        }
    }

    void print_bound(const context_ptr& ctx,
                     floatingpointT lower_bound,
                     floatingpointT upper_bound,
                     maximize_tag)
    {
        bool better_gap = (lower_bound > bestlb || upper_bound < bestub);

        if (upper_bound < bestub)
            bestub = upper_bound;

        if (lower_bound > bestlb)
            bestlb = lower_bound;

        if (better_gap) {
            if (bestlb == static_cast<floatingpointT>(0.0))
                info(ctx, "  - Upper bound: {}   (gap: 0%)\n", bestub);
            else
                info(ctx,
                     "  - Upper bound: {}   (gap: {}%)\n",
                     bestub,
                     static_cast<floatingpointT>(100.) * (bestlb - bestub) /
                       bestlb);
        }
    }

    floatingpointT init_ub(minimize_tag)
    {
        return std::numeric_limits<floatingpointT>::max();
    }

    floatingpointT init_ub(maximize_tag)
    {
        return std::numeric_limits<floatingpointT>::lowest();
    }

    template<typename SolverT>
    void operator()(const SolverT& slv,
                    const context_ptr& ctx,
941
                    const result& best)
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
    {
        floatingpointT lb = init_bound(slv);
        floatingpointT ub = init_ub(modeT());

        if (not best.solutions.empty())
            ub = static_cast<floatingpointT>(best.solutions.back().value);

        for (auto j = 0; j != slv.n; ++j)
            lb += add_bound(slv, j, slv.compute_sum_A_pi(j), modeT());

        lb *= max_cost; // restore original cost

        print_bound(ctx, lb, ub, modeT());
    }
};

958
struct compute_order
959
960
{
    std::vector<int> R;
961
962
    std::vector<std::pair<int, int>> m_order;
    solver_parameters::constraint_order order;
963
    bool use_cycle;
964

965
966
967
968
    compute_order(const solver_parameters::constraint_order order_,
                  const int variable_number)
      : R(static_cast<std::vector<int>::size_type>(variable_number))
      , m_order(static_cast<std::vector<int>::size_type>(variable_number))
969
970
971
972
      , order(order_ == solver_parameters::constraint_order::cycle
                ? solver_parameters::constraint_order::none
                : order_)
      , use_cycle(order_ == solver_parameters::constraint_order::cycle)
973
    {}
974

975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
    static solver_parameters::constraint_order next_state(
      solver_parameters::constraint_order current) noexcept
    {
        using int_type = typename std::underlying_type<
          solver_parameters::constraint_order>::type;

        auto int_current = static_cast<int_type>(current);
        auto int_max =
          static_cast<int_type>(solver_parameters::constraint_order::cycle);

        ++int_current;

        return static_cast<solver_parameters::constraint_order>(
          int_current >= int_max ? 0 : int_current);
    }

991
992
    template<typename Solver, typename Xtype>
    void init(const Solver& s, const Xtype& x)
993
    {
994
995
996
997
998
999
1000
        switch (order) {
        case solver_parameters::constraint_order::infeasibility_decr:
        case solver_parameters::constraint_order::infeasibility_incr:
            infeasibility_local_compute_violated_constraints(s, x);
            break;
        case solver_parameters::constraint_order::pi_sign_change:
            std::iota(R.begin(), R.end(), 0);
For faster browsing, not all history is shown. View entire blame