fix #567

src/algorithms/nsga2.cpp

@@ -128,6 +128,19 @@ population nsga2::evolve(population pop) const
     }
     // ---------------------------------------------------------------------------------------------------------
 
+    // We check that the initial population individuals are within the problem bounds.
+    for (decltype(NP) i = 0u; i < NP; ++i) {
+        const auto &x = pop.get_x()[i];
+        for (decltype(x.size()) j = 0u; j < x.size(); ++j) {
+            if (x[j] < bounds.first[j] || x[j] > bounds.second[j]) {
+                pagmo_throw(std::invalid_argument,
+                            "Individual " + std::to_string(i) + " has a gene at position "
+                                + std::to_string(j) + " that is outside the problem bounds. "
+                                + get_name() + " cannot deal with it.");
+            }
+        }
+    }
+
     // No throws, all valid: we clear the logs
     m_log.clear();
 

src/utils/genetic_operators.cpp

@@ -112,10 +112,8 @@ std::pair<vector_double, vector_double> sbx_crossover_impl(const vector_double &
                 betaq = sbx_betaq(beta, eta_c, rand01);
                 c2 = 0.5 * ((y1 + y2) + betaq * (y2 - y1));
 
-                if (c1 < lb[i]) c1 = lb[i];
-                if (c2 < lb[i]) c2 = lb[i];
-                if (c1 > ub[i]) c1 = ub[i];
-                if (c2 > ub[i]) c2 = ub[i];
+                c1 = std::isfinite(c1) ? std::clamp(c1, lb[i], ub[i]) : lb[i];
+                c2 = std::isfinite(c2) ? std::clamp(c2, lb[i], ub[i]) : lb[i];
                 if (drng(random_engine) < .5) {
                     child1[i] = c1;
                     child2[i] = c2;

tests/genetic_operators.cpp

@@ -120,6 +120,54 @@ BOOST_AUTO_TEST_CASE(sbx_crossover_test)
         });
 }
 
+BOOST_AUTO_TEST_CASE(sbx_crossover_nan_inf_test)
+{
+    detail::random_engine_type random_engine(32u);
+    // Parents with genes severely violating bounds can cause a negative beta in the
+    // SBX crossover calculation, which leads to NaN in std::pow(). We verify the
+    // children are always finite (no NaN/inf propagation) and clamped to bounds
+    // when crossover occurs.
+    vector_double lb = {0., 0., -5.};
+    vector_double ub = {1., 1., 5.};
+    // parent1[0] is far below lb[0], causing a negative beta in sbx_betaq
+    // which would previously produce NaN in pow().
+    vector_double parent1 = {-10., 0.3, 1.};
+    vector_double parent2 = {0.7, 0.8, -2.};
+    for (int trial = 0; trial < 100; ++trial) {
+        // Use p_cr = 1 to force crossover to always run.
+        auto children
+            = sbx_crossover(parent1, parent2, {lb, ub}, 0u, 1.0, 10., random_engine);
+        BOOST_CHECK(std::isfinite(children.first[0]));
+        BOOST_CHECK(std::isfinite(children.first[1]));
+        BOOST_CHECK(std::isfinite(children.first[2]));
+        BOOST_CHECK(std::isfinite(children.second[0]));
+        BOOST_CHECK(std::isfinite(children.second[1]));
+        BOOST_CHECK(std::isfinite(children.second[2]));
+        // Genes that were crossbred must be within bounds. A gene retains the
+        // parent value (even if out of bounds) when crossover is skipped for
+        // that gene, so we only verify that the clamped-to-bounds invariant
+        // holds when the value differs from the parent copy.
+        if (std::abs(children.first[0] - parent1[0]) > 1e-14) {
+            BOOST_CHECK(children.first[0] >= lb[0] && children.first[0] <= ub[0]);
+        }
+        if (std::abs(children.first[1] - parent1[1]) > 1e-14) {
+            BOOST_CHECK(children.first[1] >= lb[1] && children.first[1] <= ub[1]);
+        }
+        if (std::abs(children.first[2] - parent1[2]) > 1e-14) {
+            BOOST_CHECK(children.first[2] >= lb[2] && children.first[2] <= ub[2]);
+        }
+        if (std::abs(children.second[0] - parent2[0]) > 1e-14) {
+            BOOST_CHECK(children.second[0] >= lb[0] && children.second[0] <= ub[0]);
+        }
+        if (std::abs(children.second[1] - parent2[1]) > 1e-14) {
+            BOOST_CHECK(children.second[1] >= lb[1] && children.second[1] <= ub[1]);
+        }
+        if (std::abs(children.second[2] - parent2[2]) > 1e-14) {
+            BOOST_CHECK(children.second[2] >= lb[2] && children.second[2] <= ub[2]);
+        }
+    }
+}
+
 BOOST_AUTO_TEST_CASE(polynomial_mutation_test)
 {
     detail::random_engine_type random_engine(32u);

tests/nsga2.cpp

@@ -136,6 +136,39 @@ BOOST_AUTO_TEST_CASE(nsga2_evolve_test)
     pop4 = user_algo2.evolve(pop4);
 }
 
+struct mo_bounds_check {
+    vector_double fitness(const vector_double &) const
+    {
+        return {0., 0.};
+    }
+    vector_double::size_type get_nobj() const
+    {
+        return 2u;
+    }
+    std::pair<vector_double, vector_double> get_bounds() const
+    {
+        return {{0., 0.}, {1., 1.}};
+    }
+};
+
+BOOST_AUTO_TEST_CASE(nsga2_out_of_bounds_population_test)
+{
+    // Create a population with individuals outside bounds.
+    // NSGA-II should throw on evolve.
+    problem prob{mo_bounds_check{}};
+    population pop{prob, 8u, 23u};
+    // Modify the first individual to be outside the bounds.
+    auto x = pop.get_x()[0];
+    auto f = pop.get_f()[0];
+    x[0] = -1.0; // below lower bound
+    pop.set_xf(0, x, f);
+    BOOST_CHECK_THROW(nsga2{1u}.evolve(pop), std::invalid_argument);
+    // Test with a value above the upper bound.
+    x[0] = 2.0;
+    pop.set_xf(0, x, f);
+    BOOST_CHECK_THROW(nsga2{1u}.evolve(pop), std::invalid_argument);
+}
+
 BOOST_AUTO_TEST_CASE(nsga2_setters_getters_test)
 {
     nsga2 user_algo{1u, 0.95, 10., 0.01, 50., 32u};