20211252547

Committed 14 Dec 2025 05:03PM UTC coverage: 67.712% (+30.3%) from 37.371%

Build # 20211252547

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #303 from libKriging/multistart-optim_jemalloc

Multistart optim

Coverage Stats

3946 of 4387 new or added lines in 30 files covered. (89.95%)

20 existing lines in 5 files now uncovered.

8403 of 12410 relevant lines covered (67.71%)

53457.78 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

69.75

/tests/KrigingTest.cpp

// Must be before any include that might include <cmath>
#define _USE_MATH_DEFINES

#define CATCH_CONFIG_MAIN
#include "libKriging/utils/lk_armadillo.hpp"
#define CATCH_CONFIG_ENABLE_BENCHMARKING
#include <catch2/catch.hpp>
#include <chrono>
#include <cmath>
#include <fstream>
#include <sstream>
#include <vector>
#include <libKriging/Kriging.hpp>
#include <libKriging/KrigingLoader.hpp>
#include <libKriging/Trend.hpp>

auto f = [](const arma::rowvec& row) {
  double sum = 0;
  for (auto&& x : row) {
    // sum += ((x - .5) * (x - .5));  // cas 1
    sum += (x * x);  // cas 2
  }
  return sum;
};

auto prepare_and_run_bench = [](auto&& bench) {
  const int count = 11;
  const auto i = GENERATE_COPY(range(0, count));

  arma::arma_rng::seed_type seed_val = 123;  // populate somehow (fixed value => reproducible)

  const double logn = 1 + 0.1 * i;
  const arma::uword n = floor(std::pow(10., logn));
  const arma::uword d = floor(2 + i / 3);

  INFO("dimensions are n=" << n << " x d=" << d);

  arma::arma_rng::set_seed(seed_val);
  arma::mat X(n, d, arma::fill::randu);
  arma::colvec y(n);
  for (arma::uword k = 0; k < n; ++k)
    y(k) = f(X.row(k));

  bench(y, X, i);

  CHECK(true);
};

TEST_CASE("workflow") {
  prepare_and_run_bench([](const arma::colvec& y, const arma::mat& X, int) {
    Kriging ok = Kriging("gauss");
    Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
    ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", parameters);  // FIXME no move
    const double theta = 0.5;
    arma::vec theta_vec(X.n_cols);
    theta_vec.fill(theta);
    return std::get<1>(ok.logLikelihoodFun(theta_vec, true, false, false));
  });
}

TEST_CASE("save & reload") {
  prepare_and_run_bench([](const arma::colvec& y, const arma::mat& X, int) {
    Kriging ok = Kriging("gauss");
    Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
    ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", parameters);  // FIXME no move
    ok.save("dump.json");

    Kriging ok_reloaded = Kriging::load("dump.json");
    assert(KrigingLoader::describe("dump.json") == KrigingLoader::KrigingType::Kriging);

    const double theta = 0.5;
    arma::vec theta_vec(X.n_cols);
    theta_vec.fill(theta);
    return std::get<1>(ok_reloaded.logLikelihoodFun(theta_vec, true, false, false));
  });
}

TEST_CASE("fit benchmark", "[.benchmark]") {
  prepare_and_run_bench([](const arma::colvec& y, const arma::mat& X, int i) {
    Kriging ok = Kriging("gauss");
    BENCHMARK("Kriging::fit#" + std::to_string(i)) {
      Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
      return ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", parameters);  // FIXME no move
    };
  });
}

TEST_CASE("logLikelihoodFun benchmark", "[.benchmark]") {
  prepare_and_run_bench([](const arma::colvec& y, const arma::mat& X, int i) {
    Kriging ok = Kriging("gauss");
    ok.fit(y, X);  // FIXME no move

    const double theta = 0.5;
    arma::vec theta_vec(X.n_cols);
    theta_vec.fill(theta);

    BENCHMARK("Kriging::logLikelihoodFun#" + std::to_string(i)) {
      return std::get<0>(ok.logLikelihoodFun(theta_vec, false, false, false));  //
    };
  });
}

TEST_CASE("logLikelihoodGrad benchmark", "[.benchmark]") {
  prepare_and_run_bench([](const arma::colvec& y, const arma::mat& X, int i) {
    Kriging ok = Kriging("gauss");
    ok.fit(y, X);  // FIXME no move

    const double theta = 0.5;
    arma::vec theta_vec(X.n_cols);
    theta_vec.fill(theta);

    BENCHMARK("Kriging::logLikelihoodGrad#" + std::to_string(i)) {
      return std::get<1>(ok.logLikelihoodFun(theta_vec, true, false, false));  //
    };
  });
}

// Branin function for testing
auto branin = [](const arma::rowvec& x) {
  double x1 = x(0) * 15.0 - 5.0;
  double x2 = x(1) * 15.0;
  double pi = M_PI;
  double term1 = x2 - 5.0 / (4.0 * pi * pi) * (x1 * x1) + 5.0 / pi * x1 - 6.0;
  return term1 * term1 + 10.0 * (1.0 - 1.0 / (8.0 * pi)) * std::cos(x1) + 10.0;
};

TEST_CASE("Branin BFGS") {
  arma::arma_rng::set_seed(42);

  // Generate 20 uniformly distributed points in [0,1]^2
  const arma::uword n = 20;
  const arma::uword d = 2;
  arma::mat X(n, d, arma::fill::randu);

  // Compute Branin function values
  arma::colvec y(n);
  for (arma::uword k = 0; k < n; ++k) {
    y(k) = branin(X.row(k));
  }

  SECTION("BFGS single-start") {
    Kriging ok = Kriging("gauss");
    Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
    ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", parameters);

    // Check that the model was fitted successfully
    CHECK(ok.theta().n_elem == d);
    CHECK(ok.sigma2() > 0);

    // Test prediction on a new point
    arma::mat X_new(1, d);
    X_new(0, 0) = 0.5;
    X_new(0, 1) = 0.5;
    auto pred = ok.predict(X_new, true, false, false);
    CHECK(std::get<0>(pred).n_elem == 1);
    CHECK(std::get<1>(pred).n_elem == 1);
  }

  SECTION("BFGS20 multistart") {
    Kriging ok = Kriging("gauss");
    Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
    ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);

    // Check that the model was fitted successfully
    CHECK(ok.theta().n_elem == d);
    CHECK(ok.sigma2() > 0);

    // Test prediction on a new point
    arma::mat X_new(1, d);
    X_new(0, 0) = 0.5;
    X_new(0, 1) = 0.5;
    auto pred = ok.predict(X_new, true, false, false);
    CHECK(std::get<0>(pred).n_elem == 1);
    CHECK(std::get<1>(pred).n_elem == 1);
  }

  SECTION("BFGS vs BFGS20 comparison") {
    // Fit with single-start BFGS
    Kriging ok_single = Kriging("gauss");
    Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
    ok_single.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", parameters);
    double ll_single = ok_single.logLikelihood();

    // Fit with multi-start BFGS20
    Kriging ok_multi = Kriging("gauss");
    ok_multi.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);
    double ll_multi = ok_multi.logLikelihood();

    // Multi-start should find at least as good solution as single-start
    INFO("Single-start log-likelihood: " << ll_single);
    INFO("Multi-start log-likelihood: " << ll_multi);

    // Check and display OpenMP status
#ifdef _OPENMP
    INFO("OpenMP is ENABLED - parallel execution active");
#else
    INFO("OpenMP is DISABLED - sequential execution only");
    WARN("Enable OpenMP for parallel speedup: configure with -DCMAKE_CXX_FLAGS=\"-fopenmp\"");
#endif

    CHECK(ll_multi >= ll_single - std::abs(ll_single) * 1e-6);  // Allow small numerical tolerance (relative)
  }

  SECTION("BFGS10 vs BFGS20 timing (parallel efficiency)") {
    using namespace std::chrono;

    const int n_runs = 3;
    std::vector<double> times10, times20;

    INFO("Running " << n_runs << " iterations of each test...");

    // Run multiple times to get better statistics
    for (int run = 0; run < n_runs; ++run) {
      // Time BFGS10
      auto start10 = high_resolution_clock::now();
      Kriging ok10 = Kriging("gauss");
      Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
      ok10.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS10", "LL", parameters);
      auto end10 = high_resolution_clock::now();
      double time10 = duration_cast<milliseconds>(end10 - start10).count();
      times10.push_back(time10);

      // Time BFGS20
      auto start20 = high_resolution_clock::now();
      Kriging ok20 = Kriging("gauss");
      ok20.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);
      auto end20 = high_resolution_clock::now();
      double time20 = duration_cast<milliseconds>(end20 - start20).count();
      times20.push_back(time20);

      INFO("Run " << (run + 1) << ": BFGS10=" << time10 << "ms, BFGS20=" << time20 << "ms");
    }

    // Calculate averages
    double avg10 = 0.0, avg20 = 0.0;
    for (int i = 0; i < n_runs; ++i) {
      avg10 += times10[i];
      avg20 += times20[i];
    }
    avg10 /= n_runs;
    avg20 /= n_runs;

    double ratio = (avg10 > 0) ? (avg20 / avg10) : 0.0;

    INFO("Average BFGS10 time: " << avg10 << " ms");
    INFO("Average BFGS20 time: " << avg20 << " ms");
    INFO("Average ratio (time20/time10): " << ratio);

    // With effective parallelization, BFGS20 should take similar time to BFGS10
    // The ratio should be close to 1.0, certainly less than 1.5
    // Without parallelization, the ratio would be close to 2.0
    if (ratio < 1.3) {
      INFO("Ratio < 1.3: Parallelization is working effectively!");
    } else if (ratio < 1.7) {
      INFO("Ratio between 1.3-1.7: Partial parallelization or overhead");
    } else {
      INFO("Ratio >= 1.7: Likely sequential execution (no parallelization)");
    }

    // Basic sanity checks
    CHECK(avg10 > 0);
    CHECK(avg20 > 0);
    CHECK(ratio > 0);
  }

  SECTION("BFGS20 vs 20×BFGS1 equivalence (given theta0)") {
    // Test that BFGS20 with a specific theta0 matrix gives the same result
    // as running 20 separate BFGS1 optimizations, one for each row of theta0
    
    INFO("Testing that BFGS20 with given theta0 equals 20 separate BFGS1 runs");
    
    // Create a specific theta0 matrix with 20 starting points
    arma::mat theta0(20, d);
    arma::arma_rng::set_seed(12345);  // Fixed seed for reproducibility
    for (arma::uword i = 0; i < 20; i++) {
      theta0.row(i) = arma::randu<arma::rowvec>(d) % (arma::ones<arma::rowvec>(d) * 2.0) + 0.01;
    }
    
    INFO("Created theta0 matrix: " << theta0.n_rows << " x " << theta0.n_cols);
    
    // Run BFGS20 with this theta0
    INFO("Running BFGS20...");
    Kriging ok_bfgs20 = Kriging("gauss");
    Kriging::Parameters params_multi{std::nullopt, false, theta0, true, std::nullopt, true};
    ok_bfgs20.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", params_multi);
    
    arma::vec theta_bfgs20 = ok_bfgs20.theta();
    double sigma2_bfgs20 = ok_bfgs20.sigma2();
    double ll_bfgs20 = ok_bfgs20.logLikelihood();
    
    INFO("BFGS20 final result: theta=" << theta_bfgs20.t() << ", sigma2=" << sigma2_bfgs20 << ", LL=" << ll_bfgs20);
    
    // Run 20 separate BFGS1 optimizations, one for each row of theta0
    std::vector<double> ll_vec(20);
    std::vector<arma::vec> theta_vec(20);
    std::vector<double> sigma2_vec(20);
    
    INFO("Running 20 separate BFGS1 optimizations...");
    
    // Important: Set seed again before BFGS1 runs to match conditions
    arma::arma_rng::set_seed(12345);
    
    for (arma::uword i = 0; i < 20; i++) {
      Kriging ok_bfgs1 = Kriging("gauss");
      arma::mat theta0_i = theta0.row(i);  // Extract row i as 1×d matrix
      Kriging::Parameters params_single{std::nullopt, false, theta0_i, true, std::nullopt, true};
      ok_bfgs1.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", params_single);
      
      theta_vec[i] = ok_bfgs1.theta();
      sigma2_vec[i] = ok_bfgs1.sigma2();
      ll_vec[i] = ok_bfgs1.logLikelihood();
    }
    
    INFO("Completed all 20 BFGS1 runs");
    
    // Find the best result from the 20 BFGS1 runs
    auto best_iter = std::max_element(ll_vec.begin(), ll_vec.end());
    int best_idx = std::distance(ll_vec.begin(), best_iter);
    double ll_best_bfgs1 = ll_vec[best_idx];
    arma::vec theta_best_bfgs1 = theta_vec[best_idx];
    double sigma2_best_bfgs1 = sigma2_vec[best_idx];
    
    INFO("Best BFGS1[" << best_idx << "]: theta=" << theta_best_bfgs1.t() 
         << ", sigma2=" << sigma2_best_bfgs1 << ", LL=" << ll_best_bfgs1);
    
    // Compare BFGS20 result with best BFGS1 result
    // They should be very close (within numerical tolerance)
    double theta_diff = arma::norm(theta_bfgs20 - theta_best_bfgs1, 2);
    double theta_norm = arma::norm(theta_best_bfgs1, 2);
    double sigma2_diff = std::abs(sigma2_bfgs20 - sigma2_best_bfgs1);
    double ll_diff = std::abs(ll_bfgs20 - ll_best_bfgs1);
    
    INFO("Differences:");
    INFO("  ||theta_BFGS20 - theta_best_BFGS1||_2 = " << theta_diff);
    INFO("  |sigma2_BFGS20 - sigma2_best_BFGS1| = " << sigma2_diff);
    INFO("  |LL_BFGS20 - LL_best_BFGS1| = " << ll_diff);
    
    // Check for EXACT equivalence
    // BFGS20 and 20×BFGS1 should give identical results (same algorithm, same starting points)
    // With proper thread-local optimization and staggered startup, results should be deterministic
    CHECK(theta_diff / theta_norm <= 1e-3);  // Theta must be exactly the same (relative)
    CHECK(sigma2_diff / std::abs(sigma2_best_bfgs1) <= 1e-3); // Sigma2 must be exactly the same (relative)
    CHECK(ll_diff / std::abs(ll_best_bfgs1) <= 1e-3);     // Log-likelihood must be exactly the same (relative)
    
    INFO("✓ BFGS20 and best of 20×BFGS1 are EXACTLY equivalent!");
  }
}

TEST_CASE("Kriging fit with given parameters - BFGS1") {
  arma::arma_rng::set_seed(123);
  const arma::uword n = 20;
  const arma::uword d = 2;

  arma::mat X(n, d, arma::fill::randu);
  arma::colvec y(n);
  for (arma::uword k = 0; k < n; ++k)
    y(k) = f(X.row(k));

  // Define specific starting parameters
  arma::mat theta_start(1, d);
  theta_start(0, 0) = 0.5;
  theta_start(0, 1) = 0.3;
  double sigma2_start = 0.1;

  Kriging ok = Kriging("gauss");
  // Provide starting values for sigma2 and theta, optimize all
  Kriging::Parameters parameters{sigma2_start, true, theta_start, true, std::nullopt, true};
  ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", parameters);

  // Check that optimization ran (parameters should be different from starting values)
  CHECK(ok.theta().n_elem == d);
  CHECK(ok.sigma2() > 0);
  
  // Verify predictions work
  arma::mat X_new(1, d);
  X_new.fill(0.5);
  auto pred = ok.predict(X_new, true, false, false);
  CHECK(std::get<0>(pred).n_elem == 1);
  CHECK(std::get<1>(pred).n_elem == 1);
}

TEST_CASE("Kriging fit with given parameters - BFGS20") {
  arma::arma_rng::set_seed(123);
  const arma::uword n = 20;
  const arma::uword d = 2;

  arma::mat X(n, d, arma::fill::randu);
  arma::colvec y(n);
  for (arma::uword k = 0; k < n; ++k)
    y(k) = f(X.row(k));

  // Define specific theta starting points for multistart
  arma::mat theta_starts(20, d);
  for (arma::uword i = 0; i < 20; i++) {
    theta_starts.row(i) = arma::randu<arma::rowvec>(d) * 2.0 + 0.1;
  }
  double sigma2_given = 0.1;

  Kriging ok = Kriging("gauss");
  Kriging::Parameters parameters{sigma2_given, false, theta_starts, true, std::nullopt, true};
  ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);

  // Check that sigma2 was kept fixed
  CHECK(ok.sigma2() == sigma2_given);
  // Theta should be optimized (one of the starting points)
  CHECK(ok.theta().n_elem == d);

  // Verify predictions work
  arma::mat X_new(1, d);
  X_new.fill(0.5);
  auto pred = ok.predict(X_new, true, false, false);
  CHECK(std::get<0>(pred).n_elem == 1);
  CHECK(std::get<1>(pred).n_elem == 1);
}

TEST_CASE("Kriging all parameter combinations", "[multistart]") {
  arma::arma_rng::set_seed(123);
  const arma::uword n = 20;
  const arma::uword d = 2;

  arma::mat X(n, d, arma::fill::randu);
  arma::colvec y(n);
  for (arma::uword k = 0; k < n; ++k)
    y(k) = f(X.row(k));

  // Pre-defined parameter values
  double sigma2_val = 0.1;
  arma::mat theta_val(1, d);
  theta_val.fill(0.5);
  arma::mat theta_starts(20, d);
  for (arma::uword i = 0; i < 20; i++) {
    theta_starts.row(i) = arma::randu<arma::rowvec>(d) * 2.0 + 0.1;
  }

  // Test all combinations of parameter estimation flags with different optimizers
  std::vector<std::string> optims = {"none", "BFGS", "BFGS20"};
  
  for (const auto& optim : optims) {
    DYNAMIC_SECTION("Optim: " << optim) {
      // Combination 1: Estimate all parameters (not valid for "none")
      if (optim != "none") {
        SECTION("Estimate all (sigma2, theta, beta)") {
          Kriging kr("gauss");
          Kriging::Parameters params{std::nullopt, true, std::nullopt, true, std::nullopt, true};
          kr.fit(y, X, Trend::RegressionModel::Constant, false, optim, "LL", params);
          CHECK(kr.sigma2() > 0);
          CHECK(kr.theta().n_elem == d);
        }
      }

      // Combination 2: Fix sigma2, estimate theta and beta (not valid for "none")
      if (optim != "none") {
        SECTION("Fix sigma2, estimate theta and beta") {
          Kriging kr("gauss");
          Kriging::Parameters params{sigma2_val, false, std::nullopt, true, std::nullopt, true};
          kr.fit(y, X, Trend::RegressionModel::Constant, false, optim, "LL", params);
          CHECK(kr.sigma2() == sigma2_val);
          CHECK(kr.theta().n_elem == d);
        }
      }

      // Combination 3: Estimate sigma2, fix theta, estimate beta (only "none" truly fixes theta)
      SECTION("Estimate sigma2, fix theta, estimate beta") {
        Kriging kr("gauss");
        Kriging::Parameters params{std::nullopt, true, theta_val, false, std::nullopt, true};
        kr.fit(y, X, Trend::RegressionModel::Constant, false, optim, "LL", params);
        CHECK(kr.sigma2() > 0);
        if (optim == "none") {
          CHECK(arma::approx_equal(kr.theta(), theta_val.t(), "absdiff", 1e-10));
        } else {
          // BFGS variants use theta as starting point even when is_theta_estim=false
          CHECK(kr.theta().n_elem == d);
        }
      }

      // Combination 4: Fix both sigma2 and theta, estimate beta (only "none" truly fixes theta)
      SECTION("Fix sigma2 and theta, estimate beta") {
        Kriging kr("gauss");
        Kriging::Parameters params{sigma2_val, false, theta_val, false, std::nullopt, true};
        kr.fit(y, X, Trend::RegressionModel::Constant, false, optim, "LL", params);
        CHECK(kr.sigma2() == sigma2_val);
        if (optim == "none") {
          CHECK(arma::approx_equal(kr.theta(), theta_val.t(), "absdiff", 1e-10));
        } else {
          // BFGS variants use theta as starting point even when is_theta_estim=false
          CHECK(kr.theta().n_elem == d);
        }
      }

      // Combination 5: Provide starting values for sigma2 and theta (BFGS20 multistart only)
      if (optim == "BFGS20") {
        SECTION("Multistart with theta starting points") {
          Kriging kr("gauss");
          Kriging::Parameters params{sigma2_val, true, theta_starts, true, std::nullopt, true};
          kr.fit(y, X, Trend::RegressionModel::Constant, false, optim, "LL", params);
          CHECK(kr.sigma2() > 0);
          CHECK(kr.theta().n_elem == d);
        }
      }
    }
  }

  // Verify predictions work for a representative case
  Kriging kr_final("gauss");
  Kriging::Parameters params_final{std::nullopt, true, std::nullopt, true, std::nullopt, true};
  kr_final.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS", "LL", params_final);
  
  arma::mat X_new(1, d);
  X_new.fill(0.5);
  auto pred = kr_final.predict(X_new, true, false, false);
  CHECK(std::get<0>(pred).n_elem == 1);
  CHECK(std::get<1>(pred).n_elem == 1);
}

TEST_CASE("Kriging intensive stress test", "[intensive][multistart]") {
  arma::arma_rng::set_seed(789);
  
  SECTION("Many iterations with BFGS20") {
    const int n_iterations = 50;
    const arma::uword n = 30;
    const arma::uword d = 2;
    int failure_count = 0;

    INFO("Running " << n_iterations << " intensive iterations with BFGS20");
    
    for (int iter = 0; iter < n_iterations; ++iter) {
      arma::mat X = arma::randu<arma::mat>(n, d);
      arma::colvec y = arma::sin(5.0 * X.col(0)) % arma::cos(3.0 * X.col(1));
      y += 0.1 * arma::randn(n);

      Kriging ok = Kriging("gauss");
      Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
      
      try {
        ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);
        
        // Check if fit produced invalid results
        if (!std::isfinite(ok.sigma2()) || ok.sigma2() <= 0) {
          failure_count++;
          std::stringstream X_filename, y_filename;
          X_filename << "failing_kriging_X_iter" << iter << "_failure" << failure_count << ".csv";
          y_filename << "failing_kriging_y_iter" << iter << "_failure" << failure_count << ".csv";
          
          X.save(X_filename.str(), arma::csv_ascii);
          y.save(y_filename.str(), arma::csv_ascii);
          
          INFO("Saved failing case " << failure_count << " (iteration " << iter << "):");
          INFO("  X saved to: " << X_filename.str());
          INFO("  y saved to: " << y_filename.str());
          INFO("  sigma2 = " << ok.sigma2());
          
          // Still report the failure for test tracking
          CHECK(ok.sigma2() > 0);
        } else {
          CHECK(ok.theta().n_elem == d);
          CHECK(ok.sigma2() > 0);
          
          // Prediction test
          arma::mat X_new = arma::randu<arma::mat>(5, d);
          auto pred = ok.predict(X_new, true, false, false);
          CHECK(std::get<0>(pred).n_elem == 5);
        }
        
        if ((iter + 1) % 10 == 0) {
          INFO("Completed iteration " << (iter + 1) << "/" << n_iterations 
               << " (" << failure_count << " failures so far)");
        }
      } catch (const std::exception& e) {
        FAIL("Exception at iteration " << (iter + 1) << ": " << e.what());
      }
    }
    
    INFO("Total failures: " << failure_count << " out of " << n_iterations);
  }

  SECTION("Large dataset with multiple starts") {
    const arma::uword n = 100;
    const arma::uword d = 3;

    arma::mat X = arma::randu<arma::mat>(n, d);
    arma::colvec y = arma::sin(4.0 * X.col(0)) % arma::cos(3.0 * X.col(1)) % arma::exp(-2.0 * X.col(2));
    y += 0.05 * arma::randn(n);

    INFO("Testing large dataset (n=" << n << ", d=" << d << ") with BFGS20");
    
    Kriging ok = Kriging("gauss");
    Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
    
    ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);
    
    CHECK(ok.theta().n_elem == d);
    CHECK(ok.sigma2() > 0);
    
    // Large prediction batch
    arma::mat X_new = arma::randu<arma::mat>(50, d);
    auto pred = ok.predict(X_new, true, true, true);
    CHECK(std::get<0>(pred).n_elem == 50);
    CHECK(std::get<1>(pred).n_elem == 50);
    CHECK(std::get<2>(pred).n_rows == 50);
  }

  SECTION("Rapid sequential BFGS20 calls") {
    const int n_rapid = 20;
    const arma::uword n = 25;
    const arma::uword d = 2;

    INFO("Running " << n_rapid << " rapid sequential BFGS20 fits");
    
    for (int i = 0; i < n_rapid; ++i) {
      arma::mat X = arma::randu<arma::mat>(n, d);
      arma::colvec y = arma::sin(5.0 * X.col(0) + 3.0 * X.col(1));
      y += 0.05 * arma::randn(n);

      Kriging ok = Kriging("gauss");
      Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
      ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);
      
      CHECK(ok.sigma2() > 0);
      CHECK(ok.theta().n_elem == d);
    }
    
    INFO("Completed all " << n_rapid << " rapid sequential fits successfully");
  }

  SECTION("Different kernel types with BFGS20") {
    const arma::uword n = 30;
    const arma::uword d = 2;
    arma::mat X = arma::randu<arma::mat>(n, d);
    arma::colvec y = arma::sin(5.0 * X.col(0)) % arma::cos(3.0 * X.col(1));
    y += 0.1 * arma::randn(n);

    std::vector<std::string> kernels = {"gauss", "exp", "matern3_2", "matern5_2"};
    
    for (const auto& kernel : kernels) {
      INFO("Testing kernel: " << kernel);
      
      Kriging ok = Kriging(kernel);
      Kriging::Parameters parameters{std::nullopt, true, std::nullopt, true, std::nullopt, true};
      
      ok.fit(y, X, Trend::RegressionModel::Constant, false, "BFGS20", "LL", parameters);
      
      CHECK(ok.theta().n_elem == d);
      CHECK(ok.sigma2() > 0);
    }
  }
}

libKriging / libKriging / 20211252547

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous