itm-optimizer-common.hpp 27.1 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
50
51
52
53
54
55
56
57
58
59
60
61
62
/**
 * @brief Local context for each thread.
 * @details Each thread need a copy of this data to avoid mutex.
 *
 */
struct local_context
{
    random_engine rng;
    std::normal_distribution<double> choose_sol_dist;
    std::normal_distribution<double> variable_p_dist;
    std::normal_distribution<double> value_p_dist;
    std::uniform_int_distribution<int> bad_solution_choose;
    std::uniform_int_distribution<std::size_t> crossover_dist;
    std::bernoulli_distribution crossover_bastert_insertion;
63
64
65
66
67

    const double init_kappa_improve_start;
    const double init_kappa_improve_increase;
    const double init_kappa_improve_stop;

68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
    const unsigned thread_id;

    local_context(const context& ctx,
                  const unsigned thread_id_,
                  const random_engine::result_type seed)
      : rng(seed)
      , choose_sol_dist(
          ctx.parameters.init_crossover_solution_selection_mean,
          ctx.parameters.init_crossover_solution_selection_stddev)
      , variable_p_dist(ctx.parameters.init_mutation_variable_mean,
                        ctx.parameters.init_mutation_variable_stddev)
      , value_p_dist(ctx.parameters.init_mutation_value_mean,
                     ctx.parameters.init_mutation_value_stddev)
      , bad_solution_choose(ctx.parameters.init_population_size / 5,
                            ctx.parameters.init_population_size - 1)
      , crossover_dist(0)
      , crossover_bastert_insertion(
          ctx.parameters.init_crossover_bastert_insertion)
86
87
88
      , init_kappa_improve_start(ctx.parameters.init_kappa_improve_start)
      , init_kappa_improve_increase(ctx.parameters.init_kappa_improve_increase)
      , init_kappa_improve_stop(ctx.parameters.init_kappa_improve_stop)
89
90
91
92
      , thread_id(thread_id_)
    {}
};

93
template<typename Cost, typename Mode>
94
class storage
95
{
96
97
98
99
100
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>;

101
    std::vector<int> m_indices;
102
    std::vector<raw_result<Mode>> m_data;
103
    bit_array m_bastert;
104
    bit_array m_random;
105

106
107
    const Cost& costs;
    double cost_constant;
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
    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,
132
                        const std::size_t hash,
133
                        const long int loop,
134
135
                        const int remaining_constraints) noexcept
    {
136
        m_indices_reader lock{ m_indices_mutex };
137
138
139
140

        m_data[id].x = x;
        m_data[id].value = value;
        m_data[id].duration = duration;
141
        m_data[id].hash = hash;
142
143
144
        m_data[id].loop = loop;
        m_data[id].remaining_constraints = remaining_constraints;
    }
145

146
    int choose_a_bad_solution(local_context& ctx)
147
    {
148
        return ctx.bad_solution_choose(ctx.rng);
149
150
    }

151
    int choose_a_solution(local_context& ctx)
152
153
154
    {
        double value;
        do
155
            value = ctx.choose_sol_dist(ctx.rng);
156
157
158
159
160
        while (value < 0 || value > 1);

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

161
public:
162
163
164
165
166
    storage(random_engine& rng,
            const Cost& costs_,
            const double cost_constant_,
            const int population_size,
            const int variables,
167
            const std::vector<merged_constraint>& constraints_)
168
      : m_indices(population_size)
169
      , m_data(population_size)
170
      , m_bastert(variables)
171
      , m_random(variables)
172
173
      , costs(costs_)
      , cost_constant(cost_constant_)
174
      , m_size(population_size)
175
    {
176
177
        for (auto& elem : m_data)
            elem.x = bit_array(variables);
178

179
180
        init_with_bastert<Cost, Mode>(m_bastert, costs_, variables, 0);

181
        for (int i = 0, e = m_size / 2; i != e; ++i) {
182
            m_data[i].x = m_bastert;
183

184
185
            std::bernoulli_distribution dist(
              std::clamp(static_cast<double>(i) / (5 * e), 0.0, 1.0));
186

187
188
189
            for (int v = 0; v != variables; ++v)
                if (dist(rng))
                    m_data[i].x.invert(v);
190
        }
191

192
        for (int i = m_size / 2, e = m_size; i + 1 < e; i += 2) {
193
194
            init_with_random(m_data[i].x, rng, variables, 0.2);
            init_with_random(m_data[i + 1].x, rng, variables, 0.8);
195

196
197
198
199
200
201
202
203
            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));
        }
204

205
        for (int i = 0, e = m_size; i != e; ++i) {
206
            m_data[i].make_hash();
207
208
209
210
211
212
213
214
215
216
217
218
219
220
            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;
            }
221
222
        }

223
        std::iota(std::begin(m_indices), std::end(m_indices), 0);
224

225
        sort();
226
    }
227

228
    void show_population(const context& ctx) const
229
230
231
232
    {
        info(ctx, " Population {}:\n", m_indices.size());
        for (int i = 0; i != m_size; ++i)
            info(ctx,
233
                 "  - {}: value {} constraints {} hash {}\n",
234
235
                 m_indices[i],
                 m_data[m_indices[i]].value,
236
237
                 m_data[m_indices[i]].remaining_constraints,
                 m_data[m_indices[i]].hash);
238
239
    }

240
241
    void insert(local_context& ctx,
                const bit_array& x,
242
                const std::size_t hash,
243
244
                const int remaining_constraints,
                const double duration,
245
                const long int loop) noexcept
246
247
248
    {
        to_log(stdout,
               5u,
249
               "- insert advance {} (hash: {}) {}s in {} loops\n",
250
               remaining_constraints,
251
               hash,
252
253
254
               duration,
               loop);

255
        int id_to_delete = m_indices[choose_a_bad_solution(ctx)];
256

257
        to_log(stdout,
258
               5u,
259
               "- delete {} ({})\n",
260
261
               id_to_delete,
               m_data[id_to_delete].value);
262

263
        replace_result(id_to_delete,
264
265
266
                       x,
                       costs.results(x, cost_constant),
                       duration,
267
                       hash,
268
269
270
                       loop,
                       remaining_constraints);

271
        sort();
272
    }
273

274
275
    void insert(local_context& ctx,
                const bit_array& x,
276
                const std::size_t hash,
277
278
                const double value,
                const double duration,
279
                const long int loop) noexcept
280
281
    {
        to_log(stdout,
282
               5u,
283
               "- insert solution {} (hash: {}) {}s in {} loops\n",
284
               value,
285
               hash,
286
287
288
               duration,
               loop);

289
        int id_to_delete = m_indices[choose_a_bad_solution(ctx)];
290

291
        to_log(stdout,
292
               5u,
293
               "- delete {} ({})\n",
294
295
               id_to_delete,
               m_data[id_to_delete].value);
296

297
        replace_result(id_to_delete, x, value, duration, hash, loop, 0);
298

299
        sort();
300
301
    }

302
303
    bool can_be_inserted(const std::size_t hash, const int constraints) const
      noexcept
304
    {
305
        m_indices_reader lock(m_indices_mutex);
306

307
        for (int i = 0; i != m_size; ++i)
308
309
            if (m_data[i].remaining_constraints == constraints &&
                m_data[i].hash == hash)
310
311
                return false;

312
        return true;
313
    }
314

315
    bool can_be_inserted([[maybe_unused]] const std::size_t hash,
316
                         const double value) const noexcept
317
    {
318
        m_indices_reader lock(m_indices_mutex);
319

320
        for (int i = 0; i != m_size; ++i)
321
322
            if (m_data[i].remaining_constraints == 0 &&
                m_data[i].value == value && m_data[i].hash == hash)
323
                return false;
324

325
        return true;
326
    }
327

328
    const raw_result<Mode>& get_worst() const noexcept
329
    {
330
331
332
333
334
        int last = 1;
        const int end = length(m_indices);

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

336
337
        return last == end ? m_data[m_indices[last - 1]]
                           : m_data[m_indices[last]];
338
    }
339

340
341
342
343
    void get_best(int& constraints_remaining,
                  double& value,
                  double& duration,
                  long int& loop) const noexcept
344
    {
345
346
        m_indices_reader lock(m_indices_mutex);
        int id = m_indices.front();
347

348
349
350
351
        constraints_remaining = m_data[id].remaining_constraints;
        value = m_data[id].value;
        duration = m_data[id].duration;
        loop = m_data[id].loop;
352
353
    }

354
    const raw_result<Mode>& get_best(int i) const noexcept
355
    {
356
357
        return m_data[m_indices[i]];
    }
358

359
    void crossover(local_context& ctx,
360
361
362
363
364
365
366
367
368
                   bit_array& x,
                   const bit_array& first,
                   const bit_array& second)
    {
        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);

369
            x.set_block(i, ((x1 ^ x2) & ctx.crossover_dist(ctx.rng)) ^ x1);
370
371
372
        }
    }

373
    void crossover(local_context& ctx, bit_array& x)
374
    {
375
376
        m_indices_reader lock(m_indices_mutex);

377
378
        if (ctx.crossover_bastert_insertion(ctx.rng)) {
            int first = m_indices[choose_a_solution(ctx)];
379

380
            crossover(ctx, x, m_data[first].x, m_bastert);
381
382
383
384
385
386

            to_log(stdout,
                   7u,
                   "- crossover between {} ({}) and bastert\n",
                   first,
                   m_data[first].value);
387
        } else {
388
389
            int first = m_indices[choose_a_solution(ctx)];
            int second = m_indices[choose_a_solution(ctx)];
390
            while (first == second)
391
                second = m_indices[choose_a_solution(ctx)];
392

393
            crossover(ctx, x, m_data[first].x, m_data[second].x);
394
395
396
397
398
399
400
401
402

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

private:
406
    void sort() noexcept
407
408
409
410
411
    {
        m_indices_writer lock{ m_indices_mutex };

        std::sort(
          std::begin(m_indices), std::end(m_indices), [this](int i1, int i2) {
412
413
414
415
416
417
418
419
420
421
422
423
              const int cst_1 = this->m_data[i1].remaining_constraints;
              const int cst_2 = this->m_data[i2].remaining_constraints;
              const double value_1 = this->m_data[i1].value;
              const double value_2 = this->m_data[i2].value;

              if (cst_1 < cst_2)
                  return true;

              if (cst_1 == cst_2)
                  return is_better_solution<Mode>(value_1, value_2);

              return false;
424
          });
425

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

442
443
444
445
446
template<typename Cost, typename Float, typename Mode>
struct best_solution_recorder
{
    std::chrono::time_point<std::chrono::steady_clock> m_start;
    storage<Cost, Mode> m_storage;
447
    std::vector<double> m_kappa_append;
448

449
    best_solution_recorder(random_engine& rng,
450
451
452
453
454
                           const unsigned thread_number,
                           const Cost& costs,
                           const double cost_constant,
                           const int variables,
                           const std::vector<merged_constraint>& constraints,
455
                           const int population_size)
456
457
458
      : m_storage{ rng,       costs,      cost_constant, population_size,
                   variables, constraints }
      , m_kappa_append(std::size_t{ thread_number }, 0.0)
459
460
461
462
    {
        m_start = std::chrono::steady_clock::now();
    }

463
464
465
466
467
    void get_best(int& constraints_remaining,
                  double& value,
                  double& duration,
                  long int& loop) const noexcept
    {
468
        m_storage.get_best(constraints_remaining, value, duration, loop);
469
470
    }

471
    void crossover(local_context& ctx, bit_array& x)
472
    {
473
        m_storage.crossover(ctx, x);
474
    }
475

476
477
    void mutation(local_context& ctx, bit_array& x)
    {
478
        if (ctx.value_p_dist.mean() == 0.0 && ctx.value_p_dist.stddev() == 0.0)
479
480
            return;

481
482
483
        double val_p, var_p;

        do
484
            var_p = ctx.variable_p_dist(ctx.rng);
485
        while (var_p <= 0.0 || var_p >= 1.0);
486
487

        do
488
            val_p = ctx.value_p_dist(ctx.rng);
489
        while (val_p < 0.0 || val_p > 1.0);
490

491
492
493
494
495
496
        to_log(stdout,
               7u,
               "- mutation variables {}% with "
               " {}% of set\n",
               var_p,
               val_p);
497

498
499
500
        std::bernoulli_distribution dist_var_p(var_p);
        std::bernoulli_distribution dist_value_p(val_p);

501
        for (int i = 0, e = x.size(); i != e; ++i) {
502
            if (dist_var_p(ctx.rng)) {
503
504
                to_log(stdout, 9u, "- mutate variable {}\n", i);

505
                x.set(i, dist_value_p(ctx.rng));
506
507
            }
        }
508
    }
509

510
    double reinit(local_context& ctx,
511
                  const bool /*is_solution*/,
512
513
514
515
                  const double kappa_min,
                  const double kappa_max,
                  bit_array& x)
    {
516
517
518
519
        to_log(stdout, 3u, "- reinitinialization thread {}.\n", ctx.thread_id);

        double kappa = kappa_min;

520
521
        if (m_kappa_append[ctx.thread_id] < ctx.init_kappa_improve_stop) {
            m_kappa_append[ctx.thread_id] += ctx.init_kappa_improve_increase;
522
523
            kappa = kappa_min +
                    (kappa_max - kappa_min) * m_kappa_append[ctx.thread_id];
524
525
526

            to_log(stdout, 5u, "- improve with kappa {}\n", kappa);
        } else {
527
            m_kappa_append[ctx.thread_id] = ctx.init_kappa_improve_start;
528
529
530
            crossover(ctx, x);

            to_log(stdout, 5u, "- crossover\n");
531
532
        }

533
        mutation(ctx, x);
534
535

        return kappa;
536
537
    }

538
539
    void try_advance(local_context& ctx,
                     const bit_array& solution,
540
                     const int remaining_constraints,
541
                     const long int loop)
542
    {
543
544
545
        auto hash = bit_array_hash()(solution);

        if (m_storage.can_be_inserted(hash, remaining_constraints)) {
546
547
            const auto end = std::chrono::steady_clock::now();
            const auto duration = compute_duration(m_start, end);
548

549
            m_storage.insert(
550
              ctx, solution, hash, remaining_constraints, duration, loop);
551
552
553
        }
    }

554
555
    void try_update(local_context& ctx,
                    const bit_array& solution,
556
                    const double value,
557
                    const long int loop)
558
    {
559
560
561
        auto hash = bit_array_hash()(solution);

        if (m_storage.can_be_inserted(hash, value)) {
562
563
            const auto end = std::chrono::steady_clock::now();
            const auto duration = compute_duration(m_start, end);
564

565
            m_storage.insert(ctx, solution, hash, value, duration, loop);
566
567
        }
    }
568

569
    void show_population(const context& ctx) const noexcept
570
    {
571
        m_storage.show_population(ctx);
572
573
    }

574
575
576
577
578
579
    const raw_result<Mode>& get_worst() const noexcept
    {
        return m_storage.get_worst();
    }

    const raw_result<Mode>& get_best(int i) const noexcept
580
    {
581
        return m_storage.get_best(i);
582
    }
583
584
};

585
template<typename Solver, typename Float, typename Mode, typename Cost>
586
587
struct optimize_functor
{
588
589
    const context& m_ctx;
    local_context m_local_ctx;
590
591
    long int m_call_number;
    unsigned m_thread_id;
592

593
594
595
    optimize_functor(const context& ctx,
                     const unsigned thread_id,
                     const typename random_engine::result_type seed)
596
      : m_ctx{ ctx }
597
      , m_local_ctx{ ctx, thread_id, seed }
598
599
      , m_call_number{ 0 }
      , m_thread_id{ thread_id }
600
    {}
601

602
    void operator()(const std::atomic_bool& stop_task,
603
                    best_solution_recorder<Cost, Float, Mode>& best_recorder,
604
605
606
607
                    const std::vector<merged_constraint>& constraints,
                    int variables,
                    const Cost& original_costs,
                    double cost_constant)
608
    {
609
        bit_array x(variables);
610

611
        auto& p = m_ctx.parameters;
612

613
        auto norm_costs = normalize_costs<Float, Cost>(
614
          m_ctx, original_costs, m_local_ctx.rng, variables);
615

616
617
618
619
620
621
        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
622
            ? compute_delta<Float, Cost>(m_ctx, norm_costs, theta, variables)
623
624
625
626
627
628
            : 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);

629
        const auto w_limit = static_cast<long int>(p.w);
630

631
632
633
634
635
        Solver slv(m_local_ctx.rng,
                   length(constraints),
                   variables,
                   norm_costs,
                   constraints);
636

637
        compute_order compute(p.order, variables);
638
        bool is_a_solution = false;
639

640
        while (!stop_task.load()) {
641
642
            ++m_call_number;
            const auto kappa_start = static_cast<Float>(best_recorder.reinit(
643
              m_local_ctx, is_a_solution, p.kappa_min, p.kappa_max, x));
644
            auto kappa = kappa_start;
645
            compute.init(slv, x);
646

647
648
            auto best_remaining = INT_MAX;
            is_a_solution = false;
649

650
            for (long int i = 0; !stop_task.load() && i != p.limit; ++i) {
651
                auto remaining =
652
                  compute.run(slv, x, m_local_ctx.rng, kappa, delta, theta);
653
654

                if (remaining == 0) {
655
                    best_recorder.try_update(
656
657
658
659
660
                      m_local_ctx,
                      x,
                      original_costs.results(x, cost_constant),
                      i);

661
                    best_remaining = 0;
662
                    is_a_solution = true;
663
                    break;
664
665
                } else {
                    best_remaining = std::min(remaining, best_remaining);
666
                }
667

668
                if (i > w_limit)
669
670
671
672
673
674
675
676
                    kappa +=
                      kappa_step * std::pow(static_cast<Float>(remaining) /
                                              static_cast<Float>(slv.m),
                                            alpha);

                if (kappa > kappa_max)
                    break;
            }
677

678
            if (best_remaining > 0) {
679
680
                best_recorder.try_advance(
                  m_local_ctx, x, best_remaining, p.limit);
681
                continue;
682
            }
683

684
685
            for (int push = 0; !stop_task.load() && push < p.pushes_limit;
                 ++push) {
686

687
688
689
                auto remaining =
                  compute.push_and_run(slv,
                                       x,
690
                                       m_local_ctx.rng,
691
                                       pushing_k_factor,
692
693
694
695
                                       delta,
                                       theta,
                                       pushing_objective_amplifier);

696
                if (remaining == 0) {
697
                    best_recorder.try_update(
698
                      m_local_ctx,
699
700
701
                      x,
                      original_costs.results(x, cost_constant),
                      -push * p.pushing_iteration_limit - 1);
702
                    is_a_solution = true;
703
                }
704

705
                kappa = kappa_start;
706
707
                for (int iter = 0;
                     !stop_task.load() && iter < p.pushing_iteration_limit;
708
                     ++iter) {
709

710
711
                    remaining = compute.run(
                      slv, x, m_local_ctx.rng, kappa, delta, theta);
712
713

                    if (remaining == 0) {
714
                        best_recorder.try_update(
715
                          m_local_ctx,
716
                          x,
Gauthier Quesnel's avatar
Gauthier Quesnel committed
717
                          original_costs.results(x, cost_constant),
718
                          -push * p.pushing_iteration_limit - iter - 1);
719
                        is_a_solution = true;
720
721
722
723
724
725
726
727
728
729
730
                        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;
731
732
                }
            }
733
        }
734
735
736
    }
};

737
//
738
739
740
// 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.
741
742
//
inline unsigned
743
get_thread_number(const baryonyx::context& ctx) noexcept
744
745
746
{
    unsigned ret;

747
    if (ctx.parameters.thread <= 0)
748
749
        ret = std::thread::hardware_concurrency();
    else
750
        ret = static_cast<unsigned>(ctx.parameters.thread);
751
752
753
754
755
756
757

    if (ret == 0)
        return 1;

    return ret;
}

758
template<typename Solver, typename Float, typename Mode, typename Cost>
759
inline result
760
optimize_problem(const context& ctx, const problem& pb)
761
{
762
763
    result r;

764
765
    if (ctx.start)
        ctx.start(ctx.parameters);
766

767
    auto constraints{ make_merged_constraints(ctx, pb) };
768
    if (constraints.empty() || pb.vars.values.empty())
769
        return r;
770

771
    random_engine rng(init_random_generator_seed(ctx));
772

773
    auto variables = length(pb.vars.values);
774
775
    auto cost = Cost(pb.objective, variables);
    auto cost_constant = pb.objective.value;
776

777
    const auto thread = get_thread_number(ctx);
778

779
780
    std::vector<optimize_functor<Solver, Float, Mode, Cost>> functors;
    std::vector<std::thread> pool(thread);
781

782
    best_solution_recorder<Cost, Float, Mode> best_recorder(
783
784
785
786
787
788
789
      rng,
      thread,
      cost,
      cost_constant,
      variables,
      constraints,
      ctx.parameters.init_population_size);
790

791
    auto seeds = generate_seed(rng, thread);
792

793
794
    std::atomic_bool stop_task;
    stop_task.store(false);
795

796
797
    for (unsigned i = 0u; i != thread; ++i)
        functors.emplace_back(ctx, i, seeds[i]);
798

799
    for (unsigned i = 0u; i != thread; ++i)
800
        pool[i] = std::thread(std::ref(functors[i]),
801
                              std::ref(stop_task),
802
803
                              std::ref(best_recorder),
                              std::cref(constraints),
804
                              variables,
805
                              std::cref(cost),
806
                              cost_constant);
807

808
809
    const auto start = std::chrono::steady_clock::now();
    auto end = start;
810

811
812
    do {
        std::this_thread::sleep_for(std::chrono::seconds{ 1L });
813

814
        if (ctx.update) {
815
816
817
818
819
820
821
822
            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;
823

824
825
            best_recorder.get_best(
              constraints_remaining, value, duration, loop);
826

827
            ctx.update(
828
              constraints_remaining, value, loop, duration, call_number);
829
830
831
        }

        end = std::chrono::steady_clock::now();
832
    } while (!is_time_limit(ctx.parameters.time_limit, start, end));
833
834
835
836
837

    stop_task.store(true);

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

839
840
841
842
843
844
    r.strings = pb.strings;
    r.affected_vars = pb.affected_vars;
    r.variable_name = pb.vars.names;
    r.variables = variables;
    r.constraints = length(constraints);

845
    const auto& first = best_recorder.get_best(0);
846

847
848
849
    if (!first.is_solution())
        r.status = result_status::time_limit_reached;
    else
850
851
        r.status = result_status::success;

852
853
854
855
    r.duration = first.duration;
    r.loop = first.loop;
    r.remaining_constraints = first.remaining_constraints;

856
    switch (ctx.parameters.storage) {
857
858
859
860
861
862
863
864
865
    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);
866
        convert(best_recorder.get_worst(), r.solutions[1], variables);