itm-optimizer-common.hpp 28.9 KB
Newer Older
Gauthier Quesnel's avatar
2019    
Gauthier Quesnel committed
1
/* Copyright (C) 2016-2019 INRA
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 *
 * 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_OPTIMIZER_COMMON_HPP
#define ORG_VLEPROJECT_BARYONYX_SOLVER_ITM_OPTIMIZER_COMMON_HPP

#include <baryonyx/core-compare>
27
#include <baryonyx/core-utils>
28
29

#include <algorithm>
30
#include <atomic>
31
32
33
34
#include <chrono>
#include <functional>
#include <future>
#include <set>
35
#include <shared_mutex>
36
37
38
39
#include <thread>
#include <utility>
#include <vector>

40
#include "itm-common.hpp"
41
#include "private.hpp"
42
#include "result.hpp"
43
44
45
46
47
48
#include "sparse-matrix.hpp"
#include "utils.hpp"

namespace baryonyx {
namespace itm {

49
template<typename Cost, typename Mode>
50
class storage
51
{
52
53
54
55
56
57
58
59
60
private:
    mutable std::shared_mutex m_indices_mutex;
    using m_indices_writer = std::unique_lock<std::shared_mutex>;
    using m_indices_reader = std::shared_lock<std::shared_mutex>;

    mutable std::vector<std::shared_mutex> m_data_mutex;
    using m_data_writer = std::unique_lock<std::shared_mutex>;
    using m_data_reader = std::shared_lock<std::shared_mutex>;

61
    std::vector<int> m_indices;
62
    std::vector<raw_result<Mode>> m_data;
63
    bit_array m_bastert;
64
    bit_array m_random;
65

66
    std::normal_distribution<> m_choose_sol_dist;
67
68
    std::bernoulli_distribution m_crossover_bastert_insertion;

69
70
    const Cost& costs;
    double cost_constant;
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
    const int m_size;

    struct bound_indices
    {
        bound_indices(int first_, int last_)
          : first(first_)
          , last(last_)
        {}

        int first;
        int last;
    };

    bound_indices indices_bound() const noexcept
    {
        m_indices_reader lock{ m_indices_mutex };

        return bound_indices{ m_indices.front(), m_indices.back() };
    }

    void replace_result(const int id,
                        const bit_array& x,
                        const double value,
                        const double duration,
                        const std::size_t hash,
96
                        const long int loop,
97
98
99
100
101
102
103
104
105
106
107
                        const int remaining_constraints) noexcept
    {
        m_data_writer lock{ m_data_mutex[id] };

        m_data[id].x = x;
        m_data[id].value = value;
        m_data[id].duration = duration;
        m_data[id].hash = hash;
        m_data[id].loop = loop;
        m_data[id].remaining_constraints = remaining_constraints;
    }
108

109
110
111
112
113
114
115
116
117
118
    int choose_a_solution(random_engine& rng)
    {
        double value;
        do
            value = m_choose_sol_dist(rng);
        while (value < 0 || value > 1);

        return static_cast<int>(m_size * value);
    }

119
public:
120
121
122
123
124
    storage(random_engine& rng,
            const Cost& costs_,
            const double cost_constant_,
            const int population_size,
            const int variables,
125
            const std::vector<merged_constraint>& constraints_,
126
127
128
            const double crossover_bastert_insertion,
            const double crossover_solution_selection_mean,
            const double crossover_solution_selection_stddev)
129
130
131
      : m_data_mutex(population_size)
      , m_indices(population_size)
      , m_data(population_size)
132
      , m_bastert(variables)
133
      , m_random(variables)
134
135
      , m_choose_sol_dist(crossover_solution_selection_mean,
                          crossover_solution_selection_stddev)
136
      , m_crossover_bastert_insertion(crossover_bastert_insertion)
137
138
      , costs(costs_)
      , cost_constant(cost_constant_)
139
      , m_size(population_size)
140
    {
141
142
        for (auto& elem : m_data)
            elem.x = bit_array(variables);
143

144
145
        init_with_bastert<Cost, Mode>(m_bastert, costs_, variables, 0);

146
147
148
        for (int i = 0, e = m_size / 2; i != e; ++i) {
            init_with_bastert<Cost, maximize_tag>(
              m_data[i].x, costs_, variables, 0);
149

150
151
            std::bernoulli_distribution dist(
              std::clamp(static_cast<double>(i) / (5 * e), 0.0, 1.0));
152

153
154
155
            for (int v = 0; v != variables; ++v)
                if (dist(rng))
                    m_data[i].x.invert(v);
156
157
        }

158
159
160
161
162
163
164
165
166
        for (int i = m_size / 2, e = m_size; i + 1 < e; i += 2) {
            init_with_pre_solve<Cost, Mode>(
              m_data[i].x,
              m_data[i + 1].x,
              rng,
              costs_,
              constraints_,
              std::clamp(static_cast<double>(i) / (5 * e), 0.0, 1.0));
        }
167

168
        for (int i = 0, e = m_size; i != e; ++i) {
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
            m_data[i].make_hash();
            m_data[i].value = costs.results(m_data[i].x, cost_constant_);
            m_data[i].remaining_constraints = 0;

            for (int k = 0, end_k = length(constraints_); k != end_k; ++k) {
                int sum = 0;
                for (const auto& elem : constraints_[k].elements) {
                    if (m_data[i].x[elem.variable_index])
                        sum += elem.factor;
                }

                if (!(constraints_[k].min <= sum &&
                      sum <= constraints_[k].max))
                    ++m_data[i].remaining_constraints;
            }
184
185
        }

186
        std::iota(std::begin(m_indices), std::end(m_indices), 0);
187

188
        sort();
189
    }
190

191
    void insert(const bit_array& x,
192
193
194
                const std::size_t hash,
                const int remaining_constraints,
                const double duration,
195
                const long int loop) noexcept
196
    {
197
198
        auto bound = indices_bound();

199
200
        to_log(stdout,
               5u,
201
               "- insert advance {} (hash: {}) {}s in {} loops\n",
202
203
204
205
206
               remaining_constraints,
               hash,
               duration,
               loop);

207
        to_log(stdout,
208
               5u,
209
210
211
212
213
214
215
216
217
218
219
220
               "- delete {} ({})\n",
               bound.last,
               m_data[bound.last].value);

        replace_result(bound.last,
                       x,
                       costs.results(x, cost_constant),
                       duration,
                       hash,
                       loop,
                       remaining_constraints);

221
        sort();
222
    }
223

224
    void insert(const bit_array& x,
225
226
227
                const std::size_t hash,
                const double value,
                const double duration,
228
                const long int loop) noexcept
229
    {
230
231
        auto bound = indices_bound();

232
        to_log(stdout,
233
               5u,
234
               "- insert solution {} (hash: {}) {}s in {} loops\n",
235
236
237
238
239
               value,
               hash,
               duration,
               loop);

240
        to_log(stdout,
241
               5u,
242
243
244
245
246
247
               "- delete {} ({})\n",
               bound.last,
               m_data[bound.last].value);

        replace_result(bound.last, x, value, duration, hash, loop, 0);

248
        sort();
249
250
251
252
253
254
    }

    bool can_be_inserted(const std::size_t hash, const int constraints) const
      noexcept
    {
        m_indices_reader lock(m_indices_mutex);
255

256
257
258
        for (int i = 0; i != m_size; ++i) {
            const auto id = m_indices[i];
            m_data_reader lock_data(m_data_mutex[id]);
259

260
261
            if (hash == m_data[id].hash)
                return false;
262

263
264
265
266
267
            if (m_data[id].remaining_constraints > constraints)
                return true;
        }

        return false;
268
    }
269

270
271
    bool can_be_inserted(const std::size_t hash, const double value) const
      noexcept
272
    {
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
        m_indices_reader lock(m_indices_mutex);

        for (int i = 0; i != m_size; ++i) {
            const auto id = m_indices[i];
            m_data_reader lock_data(m_data_mutex[id]);

            if (hash == m_data[id].hash)
                return false;

            if (m_data[id].remaining_constraints)
                return true;

            if (is_better_solution<Mode, double>(value, m_data[id].value))
                return true;
        }

        return false;
290
    }
291

292
    const raw_result<Mode>& get_worst() const noexcept
293
    {
294
295
296
297
298
        int last = 1;
        const int end = length(m_indices);

        while (last != end && m_data[m_indices[last]].is_solution())
            ++last;
299

300
301
        return last == end ? m_data[m_indices[last - 1]]
                           : m_data[m_indices[last]];
302
    }
303

304
305
306
307
    void get_best(int& constraints_remaining,
                  double& value,
                  double& duration,
                  long int& loop) const noexcept
308
309
310
311
312
313
314
315
316
317
    {
        int id;

        {
            m_indices_reader lock(m_indices_mutex);
            id = m_indices.front();
        }

        {
            m_data_reader lock_data(m_data_mutex[id]);
318
319
320
321
            constraints_remaining = m_data[id].remaining_constraints;
            value = m_data[id].value;
            duration = m_data[id].duration;
            loop = m_data[id].loop;
322
323
324
        }
    }

325
    const raw_result<Mode>& get_best(int i) const noexcept
326
    {
327
328
        return m_data[m_indices[i]];
    }
329

330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
    void crossover(random_engine& rng,
                   bit_array& x,
                   const bit_array& first,
                   const bit_array& second)
    {
        std::uniform_int_distribution<bit_array::underlying_type> dist(0);
        std::bernoulli_distribution b_dist;

        const auto block_size = x.block_size();
        for (int i = 0; i != block_size; ++i) {
            const auto x1 = first.block(i);
            const auto x2 = second.block(i);
            const auto x_xor = x1 ^ x2;
            const auto x_rnd = dist(rng);
            const auto x_add = x_xor & x_rnd;

346
            x.set_block(i, x_add | (b_dist(rng) ? x1 : x2));
347
348
349
        }
    }

350
    void crossover(random_engine& rng, bit_array& x)
351
    {
352
        if (m_crossover_bastert_insertion(rng)) {
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
            if (m_crossover_bastert_insertion(rng)) {
                int first = m_indices[choose_a_solution(rng)];

                m_data_reader lock_data_1{ m_data_mutex[first] };
                crossover(rng, x, m_data[first].x, m_bastert);

                to_log(stdout,
                       7u,
                       "- crossover between {} ({}) and bastert\n",
                       first,
                       m_data[first].value);
            } else {
                int first = m_indices[choose_a_solution(rng)];
                init_with_random(m_random, rng, x.size(), 0.5);

                m_data_reader lock_data_1{ m_data_mutex[first] };
                crossover(rng, x, m_data[first].x, m_random);

                to_log(stdout,
                       7u,
                       "- crossover between {} ({}) and bastert\n",
                       first,
                       m_data[first].value);
            }
377
378
379
380
381
382
383
384
        } else {
            int first = m_indices[choose_a_solution(rng)];
            int second = m_indices[choose_a_solution(rng)];
            while (first == second)
                second = m_indices[choose_a_solution(rng)];

            m_data_reader lock_data_1{ m_data_mutex[first] };
            m_data_reader lock_data_2{ m_data_mutex[second] };
385
            crossover(rng, x, m_data[first].x, m_data[second].x);
386
387
388
389
390
391
392
393
394

            to_log(stdout,
                   7u,
                   "- crossover between {} ({}) and {} ({})\n",
                   first,
                   m_data[first].value,
                   second,
                   m_data[second].value);
        }
395
396
397
    }

private:
398
    void sort() noexcept
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
    {
        m_indices_writer lock{ m_indices_mutex };

        std::sort(
          std::begin(m_indices), std::end(m_indices), [this](int i1, int i2) {
              if (this->m_data[i1].remaining_constraints == 0) {
                  if (this->m_data[i2].remaining_constraints == 0)
                      return is_better_solution<Mode>(this->m_data[i1].value,
                                                      this->m_data[i2].value);
                  else
                      return true;
              } else {
                  if (this->m_data[i2].remaining_constraints == 0)
                      return false;
                  else {
                      if (this->m_data[i1].remaining_constraints ==
                          this->m_data[i2].remaining_constraints)
                          return is_better_solution<Mode>(
                            this->m_data[i1].value, this->m_data[i2].value);
                      else
                          return this->m_data[i1].remaining_constraints <
                                 this->m_data[i2].remaining_constraints;
                  }
              }
          });
424

425
426
427
428
429
430
431
432
433
434
435
436
437
#ifdef BARYONYX_ENABLE_DEBUG
        to_log(stdout, 3u, "- Solutions init population:\n");
        for (int i = 0; i != m_size; ++i) {
            to_log(stdout,
                   5u,
                   "- {} id {} value {} constraint {} hash {}\n",
                   i,
                   m_indices[i],
                   m_data[m_indices[i]].value,
                   m_data[m_indices[i]].remaining_constraints,
                   m_data[m_indices[i]].hash);
        }
#endif
438
    }
439
};
440

441
442
443
444
template<typename Cost, typename Float, typename Mode>
struct best_solution_recorder
{
    enum init_status
445
    {
446
447
448
449
450
451
452
453
454
455
456
457
458
        init_with_best = 0,
        improve_best_1,
        improve_best_2,
        improve_best_3,
        improve_best_4,
        init_with_any,
        improve_any_1,
        improve_any_2,
        improve_any_3,
        improve_any_4,
        count
    };

459
460
461
462
463
464
465
466
467
468
469
470
    static init_status advance(init_status s) noexcept
    {
        auto current = static_cast<int>(s);
        const auto count = static_cast<int>(init_status::count);

        ++current;
        if (current >= count)
            current = 0;

        return static_cast<init_status>(current);
    }

471
472
    const context_ptr& m_ctx;
    std::chrono::time_point<std::chrono::steady_clock> m_start;
473
    std::vector<init_status> m_solver_state;
474
475
476

    storage<Cost, Mode> m_storage;

477
478
    std::normal_distribution<> m_variable_p_dist;
    std::normal_distribution<> m_value_p_dist;
479
480
481
482
483
484
485

    best_solution_recorder(const context_ptr& ctx,
                           const unsigned thread_number,
                           const Cost& costs,
                           const double cost_constant,
                           const int variables,
                           const std::vector<merged_constraint>& constraints,
486
                           random_engine& rng)
487
      : m_ctx(ctx)
488
489
      , m_solver_state(thread_number, init_status::init_with_best)
      , m_storage(rng,
490
491
492
493
                  costs,
                  cost_constant,
                  ctx->parameters.init_population_size,
                  variables,
494
                  constraints,
495
496
497
498
499
500
501
                  ctx->parameters.init_crossover_bastert_insertion,
                  ctx->parameters.init_crossover_solution_selection_mean,
                  ctx->parameters.init_crossover_solution_selection_stddev)
      , m_variable_p_dist(ctx->parameters.init_mutation_variable_mean,
                          ctx->parameters.init_mutation_variable_stddev)
      , m_value_p_dist(ctx->parameters.init_mutation_value_mean,
                       ctx->parameters.init_mutation_value_stddev)
502
503
504
505
    {
        m_start = std::chrono::steady_clock::now();
    }

506
507
508
509
510
    void get_best(int& constraints_remaining,
                  double& value,
                  double& duration,
                  long int& loop) const noexcept
    {
511
        m_storage.get_best(constraints_remaining, value, duration, loop);
512
513
    }

514
515
    void mutation(random_engine& rng, bit_array& x)
    {
516
        m_storage.crossover(rng, x);
517

518
        if (m_value_p_dist.mean() == 0.0 && m_value_p_dist.stddev() == 0.0)
519
520
            return;

521
522
523
        double val_p, var_p;

        do
524
525
            var_p = m_variable_p_dist(rng);
        while (var_p <= 0.0 || var_p >= 1.0);
526
527

        do
528
529
            val_p = m_value_p_dist(rng);
        while (val_p < 0.0 || val_p > 1.0);
530

531
532
533
534
535
536
        to_log(stdout,
               7u,
               "- mutation variables {}% with "
               " {}% of set\n",
               var_p,
               val_p);
537

538
539
540
        std::bernoulli_distribution dist_var_p(var_p);
        std::bernoulli_distribution dist_value_p(val_p);

541
542
543
544
        for (int i = 0, e = x.size(); i != e; ++i) {
            if (dist_var_p(rng)) {
                to_log(stdout, 9u, "- mutate variable {}\n", i);

545
                x.set(i, dist_value_p(rng));
546
547
            }
        }
548
    }
549

550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
    double improve(const unsigned thread_id,
                   const double kappa_min,
                   const double kappa_max)
    {
        switch (m_solver_state[thread_id]) {
        case improve_best_2:
        case improve_any_2:
            return kappa_min + (kappa_max - kappa_min) * 0.1;
        case improve_best_3:
        case improve_any_3:
            return kappa_min + (kappa_max - kappa_min) * 0.15;
        case improve_best_4:
        case improve_any_4:
            return kappa_min + (kappa_max - kappa_min) * 0.20;
        default:
            return kappa_min;
        }
    }
568

569
570
    double reinit(random_engine& rng,
                  const unsigned thread_id,
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
                  const bool is_solution,
                  const double kappa_min,
                  const double kappa_max,
                  bit_array& x)
    {
        to_log(stdout,
               3u,
               "- reinitinialization thread {} - status {}. Is solution: {}\n",
               thread_id,
               m_solver_state[thread_id],
               is_solution);

        if (is_solution) {
            if (m_solver_state[thread_id] >= init_with_best &&
                m_solver_state[thread_id] < init_with_any)
                m_solver_state[thread_id] = init_with_any;
            else
                m_solver_state[thread_id] = init_with_best;
        }

        auto kappa = kappa_min;
        if (m_solver_state[thread_id] != init_with_best &&
            m_solver_state[thread_id] != init_with_any) {
            kappa = improve(thread_id, kappa_min, kappa_max);
            to_log(stdout, 5u, "- improve with kappa {}\n", kappa);
        } else {
597
            to_log(stdout, 5u, "- crossover and mutation\n");
598
            mutation(rng, x);
599
600
        }

601
        m_solver_state[thread_id] = advance(m_solver_state[thread_id]);
602
603

        return kappa;
604
605
    }

606
    void try_advance(const bit_array& solution,
607
                     const int remaining_constraints,
608
                     const long int loop)
609
    {
610
        auto hash = bit_array_hash()(solution);
611

612
613
614
        if (m_storage.can_be_inserted(hash, remaining_constraints)) {
            const auto end = std::chrono::steady_clock::now();
            const auto duration = compute_duration(m_start, end);
615

616
617
            m_storage.insert(
              solution, hash, remaining_constraints, duration, loop);
618
619
620
        }
    }

621
622
    void try_update(const bit_array& solution,
                    const double value,
623
                    const long int loop)
624
    {
625
        auto hash = bit_array_hash()(solution);
626

627
628
629
        if (m_storage.can_be_inserted(hash, value)) {
            const auto end = std::chrono::steady_clock::now();
            const auto duration = compute_duration(m_start, end);
630

631
            m_storage.insert(solution, hash, value, duration, loop);
632
633
        }
    }
634

635
636
637
638
639
640
    const raw_result<Mode>& get_worst() const noexcept
    {
        return m_storage.get_worst();
    }

    const raw_result<Mode>& get_best(int i) const noexcept
641
    {
642
        return m_storage.get_best(i);
643
    }
644
645
};

646
template<typename Solver, typename Float, typename Mode, typename Cost>
647
648
649
struct optimize_functor
{
    const context_ptr& m_ctx;
650
    random_engine m_rng;
651
652
    long int m_call_number;
    unsigned m_thread_id;
653
654
655

    optimize_functor(const context_ptr& ctx,
                     unsigned thread_id,
656
                     typename random_engine::result_type seed)
657
658
659
660
      : m_ctx{ ctx }
      , m_rng{ seed }
      , m_call_number{ 0 }
      , m_thread_id{ thread_id }
661
    {}
662

663
    void operator()(const std::atomic_bool& stop_task,
664
                    best_solution_recorder<Cost, Float, Mode>& best_recorder,
665
666
667
668
                    const std::vector<merged_constraint>& constraints,
                    int variables,
                    const Cost& original_costs,
                    double cost_constant)
669
    {
670
        bit_array x(variables);
671

672
        auto& p = m_ctx->parameters;
673

674
        auto norm_costs = normalize_costs<Float, Cost>(
675
676
          m_ctx, original_costs, m_rng, variables);

677
678
679
680
681
682
        const auto kappa_step = static_cast<Float>(p.kappa_step);
        const auto kappa_max = static_cast<Float>(p.kappa_max);
        const auto alpha = static_cast<Float>(p.alpha);
        const auto theta = static_cast<Float>(p.theta);
        const auto delta =
          p.delta < 0
Gauthier Quesnel's avatar
Gauthier Quesnel committed
683
            ? compute_delta<Float, Cost>(m_ctx, norm_costs, theta, variables)
684
685
686
687
688
689
            : static_cast<Float>(p.delta);

        const auto pushing_k_factor = static_cast<Float>(p.pushing_k_factor);
        const auto pushing_objective_amplifier =
          static_cast<Float>(p.pushing_objective_amplifier);

690
        const auto w_limit = static_cast<long int>(p.w);
691

692
        if (p.limit <= 0)
693
            p.limit = std::numeric_limits<long int>::max();
694
695
696
697
698
699

        if (p.pushes_limit <= 0)
            p.pushes_limit = 0;

        if (p.pushing_iteration_limit <= 0)
            p.pushes_limit = 0;
700

701
702
703
        Solver slv(
          m_rng, length(constraints), variables, norm_costs, constraints);

704
        compute_order compute(p.order, variables);
705
        bool is_a_solution = false;
706

707
        while (!stop_task.load()) {
708
709
710
711
            ++m_call_number;
            const auto kappa_start = static_cast<Float>(best_recorder.reinit(
              m_rng, m_thread_id, is_a_solution, p.kappa_min, p.kappa_max, x));
            auto kappa = kappa_start;
712
            compute.init(slv, x);
713

714
715
            auto best_remaining = INT_MAX;
            is_a_solution = false;
716

717
            for (long int i = 0; !stop_task.load() && i != p.limit; ++i) {
718
719
                auto remaining =
                  compute.run(slv, x, m_rng, kappa, delta, theta);
720
721

                if (remaining == 0) {
722
723
724
                    best_recorder.try_update(
                      x, original_costs.results(x, cost_constant), i);
                    best_remaining = 0;
725
                    is_a_solution = true;
726
                    break;
727
728
                } else {
                    best_remaining = std::min(remaining, best_remaining);
729
                }
730

731
                if (i > w_limit)
732
733
734
735
736
737
738
739
                    kappa +=
                      kappa_step * std::pow(static_cast<Float>(remaining) /
                                              static_cast<Float>(slv.m),
                                            alpha);

                if (kappa > kappa_max)
                    break;
            }
740

741
742
            if (best_remaining > 0) {
                best_recorder.try_advance(x, best_remaining, p.limit);
743
                continue;
744
            }
745

746
747
            for (int push = 0; !stop_task.load() && push < p.pushes_limit;
                 ++push) {
748

749
750
751
                auto remaining =
                  compute.push_and_run(slv,
                                       x,
752
                                       m_rng,
753
                                       pushing_k_factor,
754
755
756
757
                                       delta,
                                       theta,
                                       pushing_objective_amplifier);

758
                if (remaining == 0) {
759
760
761
762
                    best_recorder.try_update(
                      x,
                      original_costs.results(x, cost_constant),
                      -push * p.pushing_iteration_limit - 1);
763
                    is_a_solution = true;
764
                }
765

766
                kappa = kappa_start;
767
768
                for (int iter = 0;
                     !stop_task.load() && iter < p.pushing_iteration_limit;
769
                     ++iter) {
770

771
772
                    remaining =
                      compute.run(slv, x, m_rng, kappa, delta, theta);
773
774

                    if (remaining == 0) {
775
                        best_recorder.try_update(
776
                          x,
Gauthier Quesnel's avatar
Gauthier Quesnel committed
777
                          original_costs.results(x, cost_constant),
778
                          -push * p.pushing_iteration_limit - iter - 1);
779
                        is_a_solution = true;
780
781
782
783
784
785
786
787
788
789
790
                        break;
                    }

                    if (iter > p.w)
                        kappa +=
                          kappa_step * std::pow(static_cast<Float>(remaining) /
                                                  static_cast<Float>(slv.m),
                                                alpha);

                    if (kappa > kappa_max)
                        break;
791
792
                }
            }
793
        }
794
795
796
    }
};

797
//
798
799
800
// Get number of thread to use in optimizer from parameters list or
// from the standard thread API. If an error occurred, this function
// returns 1.
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
//
inline unsigned
get_thread_number(const baryonyx::context_ptr& ctx) noexcept
{
    unsigned ret;

    if (ctx->parameters.thread <= 0)
        ret = std::thread::hardware_concurrency();
    else
        ret = static_cast<unsigned>(ctx->parameters.thread);

    if (ret == 0)
        return 1;

    return ret;
}

818
template<typename Solver, typename Float, typename Mode, typename Cost>
819
820
821
inline result
optimize_problem(const context_ptr& ctx, const problem& pb)
{
822
823
    result r;

824
    if (ctx->start)
825
        ctx->start(ctx->parameters);
826

827
    auto constraints{ make_merged_constraints(ctx, pb) };
828
    if (constraints.empty() || pb.vars.values.empty())
829
        return r;
830

831
    random_engine rng(init_random_generator_seed(ctx));
832

833
    auto variables = length(pb.vars.values);
834
835
    auto cost = Cost(pb.objective, variables);
    auto cost_constant = pb.objective.value;
836

837
    const auto thread = get_thread_number(ctx);
838

839
840
    std::vector<optimize_functor<Solver, Float, Mode, Cost>> functors;
    std::vector<std::thread> pool(thread);
841

842
843
    best_solution_recorder<Cost, Float, Mode> best_recorder(
      ctx, thread, cost, cost_constant, variables, constraints, rng);
844

845
    auto seeds = generate_seed(rng, thread);
846

847
848
    std::atomic_bool stop_task;
    stop_task.store(false);
849

850
851
    for (unsigned i = 0u; i != thread; ++i)
        functors.emplace_back(ctx, i, seeds[i]);
852

853
    for (unsigned i = 0u; i != thread; ++i)
854
        pool[i] = std::thread(std::ref(functors[i]),
855
                              std::ref(stop_task),
856
857
                              std::ref(best_recorder),
                              std::cref(constraints),
858
                              variables,
859
                              std::cref(cost),
860
                              cost_constant);
861

862
863
    const auto start = std::chrono::steady_clock::now();
    auto end = start;
864

865
866
    do {
        std::this_thread::sleep_for(std::chrono::seconds{ 1L });
867

868
869
870
871
872
873
874
875
876
        if (ctx->update) {
            auto call_number = 0L;
            for (auto i = 0u; i != thread; ++i)
                call_number += functors[i].m_call_number;

            int constraints_remaining;
            long int loop;
            double value;
            double duration;
877

878
879
            best_recorder.get_best(
              constraints_remaining, value, duration, loop);
880

881
882
            ctx->update(
              constraints_remaining, value, loop, duration, call_number);
883
884
885
886
887
888
889
890
891
        }

        end = std::chrono::steady_clock::now();
    } while (!is_time_limit(ctx->parameters.time_limit, start, end));

    stop_task.store(true);

    for (auto& t : pool)
        t.join();
892

893
894
895
896
897
898
    r.strings = pb.strings;
    r.affected_vars = pb.affected_vars;
    r.variable_name = pb.vars.names;
    r.variables = variables;
    r.constraints = length(constraints);

899
    const auto& first = best_recorder.get_best(0);
900

901
902
903
    if (!first.is_solution())
        r.status = result_status::time_limit_reached;
    else
904
905
        r.status = result_status::success;

906
907
908
909
910
911
912
913
914
915
916
917
918
919
    r.duration = first.duration;
    r.loop = first.loop;
    r.remaining_constraints = first.remaining_constraints;

    switch (ctx->parameters.storage) {
    case solver_parameters::storage_type::one: {
        r.solutions.resize(1);
        convert(first, r.solutions[0], variables);
    } break;

    case solver_parameters::storage_type::bound: {
        r.solutions.resize(2);

        convert(first, r.solutions[0], variables);
920
        convert(best_recorder.get_worst(), r.solutions[1], variables);
921
922
923
924
925
    } break;

    case solver_parameters::storage_type::five: {
        r.solutions.resize(5);

926
927
        for (int i = 0; i != 5; ++i)
            convert(best_recorder.get_best(i), r.solutions[i], variables);
928
    } break;
929
930
    }

931
    if (ctx->finish)
932
        ctx->finish(r);
933

934
    return r;
935
936
937
938
939
940
}

} // namespace itm
} // namespace baryonyx

#endif