From 11f3c122eb2a23a4d7f8a66c4036281ff8486de0 Mon Sep 17 00:00:00 2001 From: Balaji S B Date: Tue, 18 Aug 2026 22:32:45 +0200 Subject: [PATCH 1/2] Approximations of even dimensional gaussians: closed form solutions --- .github/workflows/deploy-doxygen.yml | 1 + doxygen/doxygen_gm2d_even.h | 156 +++++++ doxygen/doxygen_main_page.h | 10 +- lib/CMakeLists.txt | 8 + .../gm_to_dirac_even_closed_form.cpp | 362 ++++++++++++++ .../gm_to_dirac_even_closed_form.h | 321 +++++++++++++ .../gm_to_dirac_even_closed_form.tpp | 132 ++++++ .../gm_to_dirac_even_optimization_params.h | 121 +++++ lib/gm_to_dirac_even/lcd_even_closed_form.h | 441 ++++++++++++++++++ plots/CMakeLists.txt | 4 + plots/generate_samples_even.cpp | 75 +++ plots/plot_samples_even.py | 253 ++++++++++ 12 files changed, 1882 insertions(+), 2 deletions(-) create mode 100644 doxygen/doxygen_gm2d_even.h create mode 100644 lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.cpp create mode 100644 lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.h create mode 100644 lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.tpp create mode 100644 lib/gm_to_dirac_even/gm_to_dirac_even_optimization_params.h create mode 100644 lib/gm_to_dirac_even/lcd_even_closed_form.h create mode 100644 plots/generate_samples_even.cpp create mode 100644 plots/plot_samples_even.py diff --git a/.github/workflows/deploy-doxygen.yml b/.github/workflows/deploy-doxygen.yml index 6310673..894e33c 100644 --- a/.github/workflows/deploy-doxygen.yml +++ b/.github/workflows/deploy-doxygen.yml @@ -37,6 +37,7 @@ jobs: GSL_RNG_SEED: "42" run: | python plots/plot_samples.py --generator build/plots/generate_samples --out-dir doxygen/images + python plots/plot_samples_even.py --generator build/plots/generate_samples_even --short-generator build/plots/generate_samples --out-dir doxygen/images - name: Build docs run: | diff --git a/doxygen/doxygen_gm2d_even.h b/doxygen/doxygen_gm2d_even.h new file mode 100644 index 0000000..543ae46 --- /dev/null +++ b/doxygen/doxygen_gm2d_even.h @@ -0,0 +1,156 @@ +/** + * \page page_gm_even gm_to_dirac_even_closed_form + * + * Gaussian-to-Dirac approximation, quadrature-free (even dimension). + * + * \section gm_even_overview Overview + * + * gm_to_dirac_even_closed_form approximates a standard normal + * distribution by a Dirac mixture with L components in N dimensions, + * exactly as gm_to_dirac_short_standard_normal_deviation does, but with + * every integral over the kernel-width variable b evaluated in closed + * form. + * + * The quadrature path issues L × N adaptive Gauss-Kronrod integrations + * per gradient evaluation. This path issues none. + * + * \section gm_even_preconditions Preconditions + * + * - N must be EVEN and >= 2. Odd N fails loudly (assert in Debug, + * approximate() returns false in Release); there is no silent fallback. + * - The target must be the standard normal N(0, I), or an isotropic + * N(0, sigma² I) through the approximate_isotropic() overloads. + * General diagonal covariance with unequal sigma_k is NOT supported: + * the scaling reduction that makes the isotropic case work does not + * exist there. + * - bMax may be ANY value > 0. Unlike the quadrature path, whose D3 uses + * a large-bMax logarithmic approximation, D3 here is the exact Ei form, + * so no "for large bMax" assumption is made anywhere. + * + * + * \section gm_even_differences Differences from the quadrature path + * + * Two behavioural differences are worth knowing when comparing results + * against gm_to_dirac_short_standard_normal_deviation: + * + * - The reported distance is the TRUE modified Cramér-von-Mises distance: + * it includes the sample-independent D1 term, which the quadrature path + * omits. + * - D3 is exact rather than logarithmically approximated. + * + * + * \section gm_even_interface Interface + * + * Inherits: + * + * - gm_to_dirac_approx_standard_normal_distribution_i + * + * Template parameter: + * + * - T ∈ {float, double}; computation is always performed in double and + * cast at the interface boundary. + * + * Provides overloads of: + * + * - approximate(...) + * - modified_van_mises_distance_sq(...) + * - modified_van_mises_distance_sq_derivative(...) + * - approximate_isotropic(...) + * - modified_van_mises_distance_sq_isotropic(...) + * + * + * \section gm_even_parameters Parameters + * + * Common parameters: + * + * - L → number of Dirac components + * - N → dimension, must be even + * - bMax → integration bound, any value > 0 + * - x → initial guess and output locations (L × N) + * - wX → weights of the Dirac mixture (optional) + * + * If wX is nullptr: + * + * - Uniform weights are assumed + * + * \section gm_even_input Supported Input Formats + * + * Three overload families are available: + * + * - Raw pointer interface (T*) + * - GSL vector interface (gsl_vector / gsl_vector_float) + * - GSL matrix interface (gsl_matrix / gsl_matrix_float) + * + * Memory layout: + * + * - x represents L × N Dirac locations, row major + * + * + * \section example_gm_even_raw Example (Raw Pointer) + * + * \code + * gm_to_dirac_even_closed_form approx; + * + * bool ok = approx.approximate( + * L, + * N, // must be even + * 100.0, // bMax, any value > 0 + * x, // initial guess / output (L × N) + * wX, // weights (optional) + * &result, + * options + * ); + * \endcode + * + * + * \section example_gm_even_distance Example (true distance) + * + * \code + * gm_to_dirac_even_closed_form approx; + * + * double distance = 0.0; + * approx.modified_van_mises_distance_sq(&distance, L, N, 100.0, x, wX); + * // distance includes D1; the quadrature path's value would be short by + * // lcd_a_n(N, 100.0) + * \endcode + * + * + * \section example_gm_even_isotropic Example (isotropic sigma² I) + * + * \code + * gm_to_dirac_even_closed_form approx; + * + * const double sigma = 2.0; + * const double bMax = 100.0; + * + * // solves the standard-normal problem at bMax / sigma and scales by sigma + * bool ok = approx.approximate_isotropic( + * L, N, sigma, bMax, x, wX, &result, options); + * \endcode + * + * + * \section gm_even_math Closed forms + * + * The distance decomposes as D = D1 - 2·D2 + D3, all expressed in units + * of pi^(N/2) — the same normalisation the quadrature path uses, so the + * two are directly comparable. The building blocks live in + * lcd_even_closed_form.h: + * + * - lcd_a_n() → D1, valid for all N (odd and even) + * - lcd_delta_bkk() → the per-sample D2 term + * - lcd_delta_bkk1() → the per-sample D2 gradient term + * - lcd_c_repulsion() → the exact per-pair D3 term + * - lcd_ei() → underflow-guarded Ei + * + * \section gm_even_notes Notes + * + * - Single-threaded implementation + * - Uses analytical gradients with no quadrature + * - The objective is assembled from reduced, O(1) terms and additionally + * carries the constant K by default + * (gm_to_dirac_even_closed_form::includeD1InObjective), so that f is the + * true distance and ApproximateOptions::ftolRel stays meaningful + * - Interface-compatible with + * gm_to_dirac_short_standard_normal_deviation + * + */ diff --git a/doxygen/doxygen_main_page.h b/doxygen/doxygen_main_page.h index bfc8734..ce1a506 100644 --- a/doxygen/doxygen_main_page.h +++ b/doxygen/doxygen_main_page.h @@ -42,8 +42,14 @@ * Implementations: * - \ref page_gm_short "gm_to_dirac_short" * - \ref page_gm_stddev "gm_to_dirac_short_standard_normal_deviation" - * - * All implementations follow the gm_to_dirac_approx_i interface. + * - \ref page_gm_even "gm_to_dirac_even_closed_form" — quadrature-free; + * requires even N and a standard normal or isotropic target, but accepts + * any bMax > 0 + * + * All implementations follow the gm_to_dirac_approx_i interface, except: + * - gm_to_dirac_short_standard_normal_deviation and + * gm_to_dirac_even_closed_form, which follow + * gm_to_dirac_approx_standard_normal_distribution_i. * * \section samples Sample Gallery * diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 15cce1e..bbe9ec3 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -36,6 +36,10 @@ set(GM_TO_DIRAC_C_WRAPPER_DIR ${GM_TO_DIRAC_DIR}/wrappers) file(GLOB GM_TO_DIRAC_SOURCES ${GM_TO_DIRAC_DIR}/*.cpp) file(GLOB GM_TO_DIRAC_C_WRAPPER_SOURCES ${GM_TO_DIRAC_C_WRAPPER_DIR}/*.cpp) +# gm_to_dirac_even sources (closed-form even-N standard normal path) +set(GM_TO_DIRAC_EVEN_DIR ${CMAKE_CURRENT_SOURCE_DIR}/gm_to_dirac_even) +file(GLOB GM_TO_DIRAC_EVEN_SOURCES ${GM_TO_DIRAC_EVEN_DIR}/*.cpp) + # common sources set(COMMON_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/debug/capture_time.cpp @@ -64,6 +68,7 @@ add_library(approxLCD_static STATIC ${DIRAC_TO_DIRAC_C_WRAPPER_SOURCES} ${GM_TO_DIRAC_SOURCES} ${GM_TO_DIRAC_C_WRAPPER_SOURCES} + ${GM_TO_DIRAC_EVEN_SOURCES} ) target_include_directories(approxLCD_static PUBLIC ${COMMON_INCLUDES} @@ -71,6 +76,7 @@ target_include_directories(approxLCD_static PUBLIC ${DIRAC_TO_DIRAC_C_WRAPPER_DIR} ${GM_TO_DIRAC_DIR} ${GM_TO_DIRAC_C_WRAPPER_DIR} + ${GM_TO_DIRAC_EVEN_DIR} ) set_target_properties(approxLCD_static PROPERTIES OUTPUT_NAME "approxLCD") target_link_libraries(approxLCD_static PUBLIC GSL::gsl GSL::gslcblas OpenMP::OpenMP_CXX) @@ -82,6 +88,7 @@ add_library(approxLCD_shared SHARED ${DIRAC_TO_DIRAC_C_WRAPPER_SOURCES} ${GM_TO_DIRAC_SOURCES} ${GM_TO_DIRAC_C_WRAPPER_SOURCES} + ${GM_TO_DIRAC_EVEN_SOURCES} ) target_include_directories(approxLCD_shared PUBLIC ${COMMON_INCLUDES} @@ -89,6 +96,7 @@ target_include_directories(approxLCD_shared PUBLIC ${DIRAC_TO_DIRAC_C_WRAPPER_DIR} ${GM_TO_DIRAC_DIR} ${GM_TO_DIRAC_C_WRAPPER_DIR} + ${GM_TO_DIRAC_EVEN_DIR} ) # platform-specific options diff --git a/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.cpp b/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.cpp new file mode 100644 index 0000000..9f71b4e --- /dev/null +++ b/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.cpp @@ -0,0 +1,362 @@ +#include "gm_to_dirac_even_closed_form.h" + +#include +#include + +#include +#include +#include + +#include "gsl_minimizer.h" +#include "gsl_utils_view_helper.h" +#include "gsl_utils_weight_helper.h" + +template +double gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + const gsl_vector* x, void* params) { + double d = 0.00; + combined_distance_metric(x, params, &d, nullptr); + return d; +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq_derivative( + const gsl_vector* x, void* params, gsl_vector* grad) { + combined_distance_metric(x, params, nullptr, grad); +} + +template +void gm_to_dirac_even_closed_form::combined_distance_metric( + const gsl_vector* x, void* params, double* f, gsl_vector* grad) { + GMToDiracEvenOptimizationParams* optiParams = + static_cast(params); + + if (f) *f = 0.00; + if (grad) gsl_vector_set_zero(grad); + + optiParams->update(x); + + calculateD2(x, optiParams, f, grad); + calculateD3(x, optiParams, f, grad); + + if (f) *f += optiParams->objectiveOffset; +} + + +template +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, double bMax, T* x, const T* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + assert(x != nullptr); + + GSLVectorView vectorViewX(x, L * N); + GSLVectorView vectorViewWX(wX, L); + return approximate(L, N, bMax, vectorViewX.get(), vectorViewWX.get(), result, + options); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + T* distance, size_t L, size_t N, double bMax, T* x, const T* wX) { + GSLVectorView vectorViewX(x, L * N); + GSLVectorView vectorViewWX(wX, L); + modified_van_mises_distance_sq(distance, L, N, bMax, vectorViewX.get(), + vectorViewWX.get()); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq_derivative( + T* gradient, size_t L, size_t N, double bMax, T* x, const T* wX) { + GSLVectorView vectorViewX(x, L * N); + GSLVectorView vectorViewWX(wX, L); + GSLVectorView vectorViewGradient(gradient, L * N); + modified_van_mises_distance_sq_derivative(vectorViewGradient.get(), L, N, + bMax, vectorViewX.get(), + vectorViewWX.get()); +} + + +template +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, double bMax, GSLMatrixType* x, const GSLVectorType* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + assert(x->size1 == L); + assert(x->size2 == N); + GSLVectorView vectorViewX(x); + return approximate(L, N, bMax, vectorViewX.get(), wX, result, options); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + T* distance, size_t L, size_t N, double bMax, GSLMatrixType* x, + const GSLVectorType* wX) { + assert(x->size1 == L); + assert(x->size2 == N); + GSLVectorView vectorViewX(x); + modified_van_mises_distance_sq(distance, L, N, bMax, vectorViewX.get(), wX); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq_derivative( + GSLMatrixType* gradient, size_t L, size_t N, double bMax, GSLMatrixType* x, + const GSLVectorType* wX) { + assert(x->size1 == L); + assert(x->size2 == N); + GSLVectorView vectorViewX(x); + GSLVectorView vectorViewGradient(gradient); + modified_van_mises_distance_sq_derivative(vectorViewGradient.get(), L, N, + bMax, vectorViewX.get(), wX); +} + + +template +bool gm_to_dirac_even_closed_form::approximate_isotropic( + size_t L, size_t N, double sigma, double bMax, T* x, const T* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + assert(sigma > 0.00); + if (!(sigma > 0.00)) return false; + + if (options.initialX) { + for (size_t i = 0; i < L * N; ++i) + x[i] = static_cast(static_cast(x[i]) / sigma); + } + + const bool success = + approximate(L, N, bMax / sigma, x, wX, result, options); + + for (size_t i = 0; i < L * N; ++i) + x[i] = static_cast(static_cast(x[i]) * sigma); + + return success; +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq_isotropic( + T* distance, size_t L, size_t N, double sigma, double bMax, T* x, + const T* wX) { + assert(sigma > 0.00); + + std::vector scaled(L * N); + for (size_t i = 0; i < L * N; ++i) + scaled[i] = static_cast(static_cast(x[i]) / sigma); + + T standardDistance = static_cast(0); + modified_van_mises_distance_sq(&standardDistance, L, N, bMax / sigma, + scaled.data(), wX); + + *distance = static_cast(sigma * sigma * + static_cast(standardDistance)); +} + + +template +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, T* x, const T* wX, GslminimizerResult* result, + const ApproximateOptions& options) { + return approximate(L, N, static_cast(options.bMax), x, wX, result, + options); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + T* distance, size_t L, size_t N, size_t bMax, T* x, const T* wX) { + modified_van_mises_distance_sq(distance, L, N, static_cast(bMax), x, + wX); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq_derivative( + T* gradient, size_t L, size_t N, size_t bMax, T* x, const T* wX) { + modified_van_mises_distance_sq_derivative(gradient, L, N, + static_cast(bMax), x, wX); +} + +template +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, GSLVectorType* x, const GSLVectorType* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + return approximate(L, N, static_cast(options.bMax), x, wX, result, + options); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + T* distance, size_t L, size_t N, size_t bMax, GSLVectorType* x, + const GSLVectorType* wX) { + modified_van_mises_distance_sq(distance, L, N, static_cast(bMax), x, + wX); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq_derivative( + GSLVectorType* gradient, size_t L, size_t N, size_t bMax, GSLVectorType* x, + const GSLVectorType* wX) { + modified_van_mises_distance_sq_derivative(gradient, L, N, + static_cast(bMax), x, wX); +} + +template +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, GSLMatrixType* x, const GSLVectorType* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + return approximate(L, N, static_cast(options.bMax), x, wX, result, + options); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + T* distance, size_t L, size_t N, size_t bMax, GSLMatrixType* x, + const GSLVectorType* wX) { + modified_van_mises_distance_sq(distance, L, N, static_cast(bMax), x, + wX); +} + +template +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq_derivative( + GSLMatrixType* gradient, size_t L, size_t N, size_t bMax, GSLMatrixType* x, + const GSLVectorType* wX) { + modified_van_mises_distance_sq_derivative(gradient, L, N, + static_cast(bMax), x, wX); +} + +// float + +template <> +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, double bMax, gsl_vector_float* x, + const gsl_vector_float* wX, GslminimizerResult* result, + const ApproximateOptions& options) { + assert(x->size == L * N); + + gsl_vector* xDouble = gsl_vector_alloc(x->size); + gsl_vector* wXDouble = nullptr; + + if (wX) { + wXDouble = gsl_vector_alloc(wX->size); + for (size_t i = 0; i < wX->size; ++i) + wXDouble->data[i] = static_cast(wX->data[i]); + } + + if (options.initialX) { + for (size_t i = 0; i < x->size; ++i) + xDouble->data[i] = static_cast(x->data[i]); + } + + gm_to_dirac_even_closed_form doubleApprox; + doubleApprox.includeD1InObjective = includeD1InObjective; + const bool success = + doubleApprox.approximate(L, N, bMax, xDouble, wXDouble, result, options); + + for (size_t i = 0; i < x->size; ++i) + x->data[i] = static_cast(xDouble->data[i]); + + gsl_vector_free(xDouble); + if (wXDouble) gsl_vector_free(wXDouble); + + return success; +} + +template <> +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + float* distance, size_t L, size_t N, double bMax, gsl_vector_float* x, + const gsl_vector_float* wX) { + double distanceDouble = 0.00; + GSLVectorView vectorViewX(x, L * N); + GSLVectorView vectorViewWX(wX, L); + + gm_to_dirac_even_closed_form doubleApprox; + doubleApprox.modified_van_mises_distance_sq(&distanceDouble, L, N, bMax, + vectorViewX.get(), + vectorViewWX.get()); + *distance = static_cast(distanceDouble); +} + +template <> +void gm_to_dirac_even_closed_form< + float>::modified_van_mises_distance_sq_derivative( + gsl_vector_float* gradient, size_t L, size_t N, double bMax, + gsl_vector_float* x, const gsl_vector_float* wX) { + gsl_vector* gradientDouble = gsl_vector_alloc(gradient->size); + + GSLVectorView vectorViewX(x, L * N); + GSLVectorView vectorViewWX(wX, L); + + gm_to_dirac_even_closed_form doubleApprox; + doubleApprox.modified_van_mises_distance_sq_derivative( + gradientDouble, L, N, bMax, vectorViewX.get(), vectorViewWX.get()); + + for (size_t i = 0; i < gradient->size; ++i) + gradient->data[i] = static_cast(gradientDouble->data[i]); + + gsl_vector_free(gradientDouble); +} + +// double + +template <> +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, double bMax, gsl_vector* x, const gsl_vector* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + assert(x != nullptr); + assert(x->size == L * N); + if (!preconditionsHold(L, N, bMax)) return false; + + if (!options.initialX) { + gsl_rng_env_setup(); + gsl_rng* r = gsl_rng_alloc(gsl_rng_default); + for (size_t i = 0; i < L; ++i) { + for (size_t d = 0; d < N; ++d) + x->data[i * N + d] = gsl_ran_gaussian(r, 1.00); // standard normal + } + gsl_rng_free(r); + } + + GSLWeightHelper wXHelper(wX, L); + GMToDiracEvenOptimizationParams params(wXHelper.get(), N, L, bMax, + includeD1InObjective); + + gsl_minimizer gslMinimizer( + options.maxIterations, options.xtolAbs, options.xtolRel, options.ftolAbs, + options.ftolRel, options.gtol, ¶ms, modified_van_mises_distance_sq, + modified_van_mises_distance_sq_derivative, combined_distance_metric); + + const int status = gslMinimizer.minimize(x, result, options.verbose); + + correctMean(x, params.wX, L, N); + + return status == GSL_SUCCESS; +} + +template <> +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + double* distance, size_t L, size_t N, double bMax, gsl_vector* x, + const gsl_vector* wX) { + assert(distance != nullptr); + if (!preconditionsHold(L, N, bMax)) return; + + GSLWeightHelper wXHelper(wX, L); + // always report the TRUE distance, offset included, regardless of what the + // minimizer's objective was configured to carry + GMToDiracEvenOptimizationParams optiParams(wXHelper.get(), N, L, bMax, true); + + *distance = modified_van_mises_distance_sq(x, &optiParams); +} + +template <> +void gm_to_dirac_even_closed_form< + double>::modified_van_mises_distance_sq_derivative(gsl_vector* gradient, + size_t L, size_t N, + double bMax, + gsl_vector* x, + const gsl_vector* wX) { + assert(gradient != nullptr); + if (!preconditionsHold(L, N, bMax)) return; + + GSLWeightHelper wXHelper(wX, L); + GMToDiracEvenOptimizationParams optiParams(wXHelper.get(), N, L, bMax, true); + + modified_van_mises_distance_sq_derivative(x, &optiParams, gradient); +} + +template class gm_to_dirac_even_closed_form; +template class gm_to_dirac_even_closed_form; diff --git a/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.h b/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.h new file mode 100644 index 0000000..2167d6e --- /dev/null +++ b/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.h @@ -0,0 +1,321 @@ +#ifndef GM_TO_DIRAC_EVEN_CLOSED_FORM_H +#define GM_TO_DIRAC_EVEN_CLOSED_FORM_H + +#include "gm_to_dirac_approx_standard_normal_distribution_i.h" +#include "gm_to_dirac_even_optimization_params.h" + +/** + * @brief quadrature-free Gaussian-to-Dirac approximation for an + * even-dimensional standard normal target + * + * Same objective as gm_to_dirac_short_standard_normal_deviation, but every + * integral over the kernel-width variable b is evaluated in closed form (see + * lcd_even_closed_form.h). + * + * One further differences from the quadrature path: + * + * - the reported distance is the TRUE mCvM distance: it includes the D1 term, + * which gm_to_dirac_short omits. + * + * @note PRECONDITIONS: the target must be the STANDARD normal N(0, I) (or an + * isotropic sigma^2 * I via approximate_isotropic()), and N must be EVEN and + * >= 2. Odd N fails loudly rather than falling back. + * + */ +template +class gm_to_dirac_even_closed_form + : public gm_to_dirac_approx_standard_normal_distribution_i { + public: + using GSLVectorType = typename gm_to_dirac_approx_standard_normal_distribution_i< + T>::GSLVectorType; + using GSLMatrixType = typename gm_to_dirac_approx_standard_normal_distribution_i< + T>::GSLMatrixType; + + /** + * @brief carry the sample-independent constant inside the minimizer's + * objective + * + * The constant does not depend on x, so omitting it leaves the gradient and + * the minimum unchanged and saves one scalar evaluation per objective call. + * With it, the objective IS the distance. + */ + bool includeD1InObjective = true; + + // clang-format off + /** + * @brief approximate a standard normal by L Dirac components + * + * @param L number of Dirac components + * @param N dimension, must be even and >= 2 + * @param bMax upper bound of the b-integral, must be > 0 + * @param x initial guess and output locations (L * N) + * @param wX weights, nullptr for uniform + * @param result minimizer result + * @param options minimizer options; options.bMax is ignored here + * @return true on success, false otherwise (including odd N) + */ + bool approximate(size_t L, + size_t N, + double bMax, + T* x, + const T* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}); + + /** + * @brief true mCvM distance, D1 included + * + * @param distance output distance + * @param L number of Dirac components + * @param N dimension, must be even and >= 2 + * @param bMax upper bound of the b-integral, must be > 0 + * @param x sample locations (L * N) + * @param wX weights, nullptr for uniform + */ + void modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, + double bMax, + T* x, + const T* wX); + + /** + * @brief gradient of the true mCvM distance with respect to x + * + * @param gradient output gradient (L * N) + * @param L number of Dirac components + * @param N dimension, must be even and >= 2 + * @param bMax upper bound of the b-integral, must be > 0 + * @param x sample locations (L * N) + * @param wX weights, nullptr for uniform + */ + void modified_van_mises_distance_sq_derivative(T* gradient, + size_t L, + size_t N, + double bMax, + T* x, + const T* wX); + + bool approximate(size_t L, + size_t N, + double bMax, + GSLVectorType* x, + const GSLVectorType* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}); + + void modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, + double bMax, + GSLVectorType* x, + const GSLVectorType* wX); + + void modified_van_mises_distance_sq_derivative(GSLVectorType* gradient, + size_t L, + size_t N, + double bMax, + GSLVectorType* x, + const GSLVectorType* wX); + + bool approximate(size_t L, + size_t N, + double bMax, + GSLMatrixType* x, + const GSLVectorType* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}); + + void modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, + double bMax, + GSLMatrixType* x, + const GSLVectorType* wX); + + void modified_van_mises_distance_sq_derivative(GSLMatrixType* gradient, + size_t L, + size_t N, + double bMax, + GSLMatrixType* x, + const GSLVectorType* wX); + // clang-format on + + // isotropic sigma^2 * I (scaling reduction) + + // clang-format off + /** + * @brief approximate an isotropic Gaussian N(0, sigma^2 * I) + * + * Solves the standard-normal problem at bMax / sigma and scales the result + * by sigma. This is why the new API takes a double bMax: bMax / sigma is + * generally not an integer. + * + * @param L number of Dirac components + * @param N dimension, must be even and >= 2 + * @param sigma standard deviation, must be > 0 + * @param bMax upper bound of the b-integral, must be > 0 + * @param x initial guess and output locations (L * N) + * @param wX weights, nullptr for uniform + * @param result minimizer result + * @param options minimizer options; options.bMax is ignored here + * @return true on success, false otherwise + * @note only isotropic covariance is supported; general diagonal covariance + * admits no such reduction + */ + bool approximate_isotropic(size_t L, + size_t N, + double sigma, + double bMax, + T* x, + const T* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}); + + /** + * @brief true mCvM distance for an isotropic Gaussian N(0, sigma^2 * I) + * + * D(sigma, bMax) = sigma^2 * D(1, bMax / sigma). + * + * @param distance output distance + * @param L number of Dirac components + * @param N dimension, must be even and >= 2 + * @param sigma standard deviation, must be > 0 + * @param bMax upper bound of the b-integral, must be > 0 + * @param x sample locations (L * N) + * @param wX weights, nullptr for uniform + */ + void modified_van_mises_distance_sq_isotropic(T* distance, + size_t L, + size_t N, + double sigma, + double bMax, + T* x, + const T* wX); + + bool approximate(size_t L, + size_t N, + T* x, + const T* wX, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}) override; + + void modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, + size_t bMax, + T* x, + const T* wX) override; + + void modified_van_mises_distance_sq_derivative(T* gradient, + size_t L, + size_t N, + size_t bMax, + T* x, + const T* wX) override; + + bool approximate(size_t L, + size_t N, + GSLVectorType* x, + const GSLVectorType* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}) override; + + void modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, + size_t bMax, + GSLVectorType* x, + const GSLVectorType* wX) override; + + void modified_van_mises_distance_sq_derivative(GSLVectorType* gradient, + size_t L, + size_t N, + size_t bMax, + GSLVectorType* x, + const GSLVectorType* wX) override; + + bool approximate(size_t L, + size_t N, + GSLMatrixType* x, + const GSLVectorType* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}) override; + + void modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, + size_t bMax, + GSLMatrixType* x, + const GSLVectorType* wX) override; + + void modified_van_mises_distance_sq_derivative(GSLMatrixType* gradient, + size_t L, + size_t N, + size_t bMax, + GSLMatrixType* x, + const GSLVectorType* wX) override; + // clang-format on + + private: + static double modified_van_mises_distance_sq(const gsl_vector* x, + void* params); + static void modified_van_mises_distance_sq_derivative(const gsl_vector* x, + void* params, + gsl_vector* grad); + static void combined_distance_metric(const gsl_vector* x, void* params, + double* f, gsl_vector* grad); + + static inline void calculateD2(const gsl_vector* x, + GMToDiracEvenOptimizationParams* params, + double* f, gsl_vector* grad); + + static inline void calculateD3(const gsl_vector* x, + GMToDiracEvenOptimizationParams* params, + double* f, gsl_vector* grad); + + static inline void correctMean(gsl_vector* x, const gsl_vector* wX, size_t L, + size_t N); + + static inline bool preconditionsHold(size_t L, size_t N, double bMax); +}; + +#include "gm_to_dirac_even_closed_form.tpp" + +template <> +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, double bMax, gsl_vector_float* x, + const gsl_vector_float* wX, GslminimizerResult* result, + const ApproximateOptions& options); + +template <> +bool gm_to_dirac_even_closed_form::approximate( + size_t L, size_t N, double bMax, gsl_vector* x, const gsl_vector* wX, + GslminimizerResult* result, const ApproximateOptions& options); + +template <> +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + float* distance, size_t L, size_t N, double bMax, gsl_vector_float* x, + const gsl_vector_float* wX); + +template <> +void gm_to_dirac_even_closed_form::modified_van_mises_distance_sq( + double* distance, size_t L, size_t N, double bMax, gsl_vector* x, + const gsl_vector* wX); + +template <> +void gm_to_dirac_even_closed_form< + float>::modified_van_mises_distance_sq_derivative( + gsl_vector_float* gradient, size_t L, size_t N, double bMax, + gsl_vector_float* x, const gsl_vector_float* wX); + +template <> +void gm_to_dirac_even_closed_form< + double>::modified_van_mises_distance_sq_derivative( + gsl_vector* gradient, size_t L, size_t N, double bMax, gsl_vector* x, + const gsl_vector* wX); + +extern template class gm_to_dirac_even_closed_form; +extern template class gm_to_dirac_even_closed_form; + +#endif // GM_TO_DIRAC_EVEN_CLOSED_FORM_H diff --git a/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.tpp b/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.tpp new file mode 100644 index 0000000..7076eff --- /dev/null +++ b/lib/gm_to_dirac_even/gm_to_dirac_even_closed_form.tpp @@ -0,0 +1,132 @@ +#ifndef GM_TO_DIRAC_EVEN_CLOSED_FORM_TPP +#define GM_TO_DIRAC_EVEN_CLOSED_FORM_TPP + +/** + * @brief check the preconditions of the even-N closed-form path + * + * @param L number of Dirac components + * @param N dimension + * @param bMax upper bound of the b-integral + * @return true if the closed form applies + */ +template +inline bool gm_to_dirac_even_closed_form::preconditionsHold(size_t L, + size_t N, + double bMax) { + assert(N % 2 == 0); + assert(N >= 2); + assert(bMax > 0.00); + assert(L >= 1); + + return N >= 2 && N % 2 == 0 && bMax > 0.00 && L >= 1; +} + +/** + * @brief attraction term: contributes -2 * D2 to f and its gradient to grad + * + * f += -2 * pow(2,k) * sum_i w_i * dBkk(k, bMax, c_i) + * + * grad += pow(2,k+1) * w_q * s_q * dBkk1(k, bMax, c_q) + * + * @param x sample locations + * @param params optimization parameters, caches must be up to date + * @param f objective accumulator, may be nullptr + * @param grad gradient accumulator, may be nullptr + */ +template +inline void gm_to_dirac_even_closed_form::calculateD2( + const gsl_vector* x, GMToDiracEvenOptimizationParams* params, double* f, + gsl_vector* grad) { + const size_t L = params->L; + const size_t N = params->N; + const size_t k = params->k; + const double bMax = params->bMax; + const gsl_vector* wX = params->wX; + + if (f) { + double sum = 0.00; + for (size_t i = 0; i < L; ++i) + sum += wX->data[i] * lcd_delta_bkk_reduced(k, bMax, params->cSqrdNorm[i]); + + *f += -2.00 * std::ldexp(sum, static_cast(k)); + } + + if (grad) { + for (size_t q = 0; q < L; ++q) { + const double factor = + std::ldexp(wX->data[q] * lcd_delta_bkk1(k, bMax, params->cSqrdNorm[q]), + static_cast(k) + 1); + for (size_t d = 0; d < N; ++d) + grad->data[q * N + d] += factor * x->data[q * N + d]; + } + } +} + +/** + * @brief repulsion term: contributes D3 to f and its gradient to grad + * + * f += sum_{i,j} w_i * w_j * (C(bMax, T_ij) - bMax^2/2) + * + * grad += 0.5 * w_q * sum_{i != q} w_i * (s_q - s_i) * Ei(-T_qi/(4 bMax^2)) + * + * @param x sample locations + * @param params optimization parameters + * @param f objective accumulator, may be nullptr + * @param grad gradient accumulator, may be nullptr + */ +template +inline void gm_to_dirac_even_closed_form::calculateD3( + const gsl_vector* x, GMToDiracEvenOptimizationParams* params, double* f, + gsl_vector* grad) { + const size_t L = params->L; + const size_t N = params->N; + const double bMax = params->bMax; + const gsl_vector* wX = params->wX; + + const double fourBMaxSqrd = 4.00 * bMax * bMax; + double d3 = 0.00; + + for (size_t i = 0; i < L; ++i) { + const double wXi = wX->data[i]; + + for (size_t j = 0; j < L; ++j) { + const double wXiwXj = wXi * wX->data[j]; + const double localDistSq = params->distanceSq(i, j); + + // coincident samples + if (localDistSq <= 0.00) continue; + + const double z = -localDistSq / fourBMaxSqrd; + const double ei = lcd_ei(z); + + if (f) + d3 += wXiwXj * (0.50 * bMax * bMax * std::expm1(z) + + 0.125 * localDistSq * ei); + + if (!grad) continue; + + const double constFactor = 0.50 * wXiwXj * ei; + for (size_t d = 0; d < N; ++d) + grad->data[i * N + d] += + constFactor * (x->data[i * N + d] - x->data[j * N + d]); + } + } + + if (f) *f += d3; +} + +template +inline void gm_to_dirac_even_closed_form::correctMean(gsl_vector* x, + const gsl_vector* wX, + size_t L, size_t N) { + std::vector mean(N, 0.00); + for (size_t i = 0; i < L; ++i) { + const double wXi = wX->data[i]; + for (size_t d = 0; d < N; ++d) mean[d] += wXi * x->data[i * N + d]; + } + for (size_t i = 0; i < L; ++i) { + for (size_t d = 0; d < N; ++d) x->data[i * N + d] -= mean[d]; + } +} + +#endif // GM_TO_DIRAC_EVEN_CLOSED_FORM_TPP diff --git a/lib/gm_to_dirac_even/gm_to_dirac_even_optimization_params.h b/lib/gm_to_dirac_even/gm_to_dirac_even_optimization_params.h new file mode 100644 index 0000000..dd6a1a5 --- /dev/null +++ b/lib/gm_to_dirac_even/gm_to_dirac_even_optimization_params.h @@ -0,0 +1,121 @@ +#ifndef GM_TO_DIRAC_EVEN_OPTIMIZATION_PARAMS_H +#define GM_TO_DIRAC_EVEN_OPTIMIZATION_PARAMS_H + +#include +#include + +#include +#include + +#include "gsl_minimizer.h" +#include "lcd_even_closed_form.h" +#include "squared_euclidean_distance_utils.h" + +/** + * @brief optimization parameters for the closed-form even-N Gaussian-to-Dirac + * approximation + */ +struct GMToDiracEvenOptimizationParams : public GslMinimizerOptimizationParams { + public: + /** + * @brief Construct a new GMToDiracEvenOptimizationParams object + * + * @param wX weights of the Dirac mixture, must be non-null and of size L + * @param N dimension of the data, must be even and >= 2 + * @param L number of Dirac components, must be >= 1 + * @param bMax upper bound of the b-integral, must be > 0 + * @param includeD1 true to carry the sample-independent constant in the + * objective (see gm_to_dirac_even_closed_form::includeD1InObjective) + */ + GMToDiracEvenOptimizationParams(const gsl_vector* wX, size_t N, size_t L, + double bMax, bool includeD1) + : GslMinimizerOptimizationParams(L, N), + wX(wX), + bMax(bMax), + k(N / 2), + reportedOffset( + lcd_reported_offset(N, bMax, weightSumOf(wX, L))), + objectiveOffset(includeD1 ? reportedOffset : 0.00), + cSqrdNorm(L, 0.00) { + assert(wX != nullptr); + assert(wX->size == L); + assert(N % 2 == 0); + assert(N >= 2); + assert(bMax > 0.00); + assert(L >= 1); + + if (L > 1) + squaredEuclideanDistanceUtilLL = + new SquaredEuclideanDistance_LL_vectorized(L, N); + } + + /** + * @brief Destroy the GMToDiracEvenOptimizationParams object + */ + ~GMToDiracEvenOptimizationParams() { + if (squaredEuclideanDistanceUtilLL) delete squaredEuclideanDistanceUtilLL; + } + + /** + * @brief refresh the per-evaluation caches for a new x + * + * Must be called whenever x changes, before reading cSqrdNorm or + * distanceSq(). + * + * @param x sample locations, L * N + */ + inline void update(const gsl_vector* x) { + assert(x->size == L * N); + + for (size_t i = 0; i < L; ++i) { + double sum = 0.00; + for (size_t d = 0; d < N; ++d) { + const double xid = x->data[i * N + d]; + sum += xid * xid; + } + cSqrdNorm[i] = sum; + } + + if (!squaredEuclideanDistanceUtilLL) return; + + const gsl_matrix xMatrix = gsl_matrix_view_array(x->data, L, N).matrix; + squaredEuclideanDistanceUtilLL->calculateDistance(&xMatrix, nullptr); + } + + /** + * @brief cached squared distance T_ij between samples i and j + * + * @param i index of the first sample + * @param j index of the second sample + * @return ||s_i - s_j||^2 + */ + inline double distanceSq(size_t i, size_t j) const { + if (!squaredEuclideanDistanceUtilLL) return 0.00; // L == 1 + return squaredEuclideanDistanceUtilLL->getDistance(i, j); + } + + const gsl_vector* wX; + const double bMax; + const size_t k; ///< N / 2 + + /// constant turning the reduced objective into the true distance: + /// distance = f_reduced + reportedOffset. See lcd_reported_offset(). + const double reportedOffset; + + /// what the minimizer's objective actually carries: reportedOffset, or 0 + /// if the caller opted out. Adding a constant costs no x-dependent + /// precision. + const double objectiveOffset; + + std::vector cSqrdNorm; ///< c_i = ||s_i||^2, size L + SquaredEuclideanDistanceUtilsLL* squaredEuclideanDistanceUtilLL = nullptr; + + private: + static double weightSumOf(const gsl_vector* wX, size_t L) { + double sum = 0.00; + for (size_t i = 0; i < L; ++i) sum += wX->data[i]; + return sum; + } +}; + +#endif // GM_TO_DIRAC_EVEN_OPTIMIZATION_PARAMS_H diff --git a/lib/gm_to_dirac_even/lcd_even_closed_form.h b/lib/gm_to_dirac_even/lcd_even_closed_form.h new file mode 100644 index 0000000..8d161b3 --- /dev/null +++ b/lib/gm_to_dirac_even/lcd_even_closed_form.h @@ -0,0 +1,441 @@ +#ifndef LCD_EVEN_CLOSED_FORM_H +#define LCD_EVEN_CLOSED_FORM_H + +#include +#include + +#include +#include +#include + +#include "gsl_quadrature_adaptive_gauss_kronrod.h" + +/** + * @file lcd_even_closed_form.h + * @brief quadrature-free closed forms for the LCD / modified Cramer-von-Mises + * distance of an even-dimensional standard normal target + * + * The mCvM distance of a Dirac mixture with samples s_i and weights w_i + * against N(0, I) decomposes into + * + * D = D1 - 2 * D2 + D3 + * + * where D1 is sample-independent, D2 is the attraction term (samples towards + * the mean) and D3 is the repulsion term (samples away from each other). + * + * @note Preconditions for the D2 helpers (lcd_delta_bkk, lcd_delta_bkk1): + * the target must be the STANDARD normal N(0, I) and N = 2k must be EVEN. + * lcd_a_n and lcd_c_repulsion are valid for all N. + */ + +/** + * @brief exponential integral Ei(x) for negative arguments + * + * @param x argument, must be strictly negative + * @return Ei(x), or exactly 0.0 in case of underflow, domain error or overflow + */ +inline double lcd_ei(double x) { + assert(x == 0.00); + + gsl_sf_result result; + const int status = gsl_sf_expint_Ei_e(x, &result); + + if (status == GSL_EUNDRFLW || status == GSL_EDOM || status == GSL_EOVRFLW) return 0.00; + assert(status == GSL_SUCCESS); + + return result.val; +} + +/** + * @brief D1, the sample-independent term of the mCvM distance + * + * Evaluates A_N(b) = integral_0^b x^(N+1) / (1 + x^2)^(N/2) dx via the + * recursion + * + * A_1(b) = 0.5 * (b * sqrt(1 + b^2) - asinh(b)) + * + * A_2(b) = 0.5 * (b^2 - log(1 + b^2)) + * + * A_n(b) = (n * A_{n-2}(b) - b^n / (1 + b^2)^((n-2)/2)) / (n - 2) + * + * applied iteratively. Valid for all N, odd and even. + * + * @param N dimension, must be >= 1 + * @param b upper integration bound, must be >= 0 + * @return A_N(b) + */ +inline double lcd_a_n(size_t N, double b) { + assert(N >= 1); + assert(b >= 0.00); + + const double bSqrd = b * b; + + size_t n = (N % 2 == 0) ? 2 : 1; + double aN = (N % 2 == 0) + ? 0.50 * (bSqrd - std::log1p(bSqrd)) + : 0.50 * (b * std::sqrt(1.00 + bSqrd) - std::asinh(b)); + + for (n += 2; n <= N; n += 2) { + const double dn = static_cast(n); + const double tail = + bSqrd * std::pow(bSqrd / (1.00 + bSqrd), (dn - 2.00) / 2.00); + aN = (dn * aN - tail) / (dn - 2.00); + } + + return aN; +} + +/** + * @brief C(b, c), the per-pair repulsion term of D3 + * + * C(b, c) = integral_0^b x * exp(-c / (4 x^2)) dx + * + * = (b^2 / 2) * exp(-c / (4 b^2)) + (c / 8) * Ei(-c / (4 b^2)) + * + * @param b upper integration bound, must be > 0 + * @param c squared distance between the two samples, must be >= 0 + * @return C(b, c); C(b, 0) = b^2 / 2 + */ +inline double lcd_c_repulsion(double b, double c) { + assert(b > 0.00); + assert(c >= 0.00); + + const double bSqrd = b * b; + if (c <= 0.00) return 0.50 * bSqrd; // coincident samples + + const double z = -c / (4.00 * bSqrd); + return 0.50 * bSqrd * std::exp(z) + 0.125 * c * lcd_ei(z); +} + +/** + * @brief C(b, c) - b^2/2, the repulsion term with its leading term removed + * + * Removed from D3 in total: (sum_i w_i)^2 * bMax^2 / 2. + * + * @param b upper integration bound, must be > 0 + * @param c squared distance between the two samples, must be >= 0 + * @return C(b, c) - b^2/2; exactly 0 at c = 0, so coincident samples + * (including every i == j pair) now contribute nothing + */ +inline double lcd_c_repulsion_reduced(double b, double c) { + assert(b > 0.00); + assert(c >= 0.00); + + if (c <= 0.00) return 0.00; // coincident samples + + const double bSqrd = b * b; + const double z = -c / (4.00 * bSqrd); + return 0.50 * bSqrd * std::expm1(z) + 0.125 * c * lcd_ei(z); +} + +/** + * @brief parameters of the raw b-integrand used by the small-c fallback + */ +struct LcdDeltaBIntegrandParams { + size_t k; ///< half the dimension + size_t q; ///< exponent of (1 + 2b^2); k for dBkk, k+1 for dBkk1 + double cValue; ///< squared norm of the sample + bool reduced; ///< true to use expm1, giving dBkk(c) - dBkk(0) directly +}; + +/** + * @brief b^(2k+1) * (1 + 2b^2)^(-q) * exp(-c / (2 (1 + 2b^2))) + * + * With reduced set, exp is replaced by expm1, so the integral evaluates + * dBkk(c) - dBkk(0) directly and never forms the O(bMax^2) intermediate - + * the same trick gm_to_dirac_short::calculateP2 uses. + * + * @param b integration variable + * @param params LcdDeltaBIntegrandParams + * @return value of the integrand at b + */ +inline double lcd_delta_b_integrand(double b, void* params) { + const LcdDeltaBIntegrandParams* p = + static_cast(params); + + const double u = 1.00 + 2.00 * b * b; + double value = b * std::pow(b * b / u, static_cast(p->k)); + for (size_t i = p->k; i < p->q; ++i) value /= u; + + const double exponent = -p->cValue / (2.00 * u); + return value * (p->reduced ? std::expm1(exponent) : std::exp(exponent)); +} + +/** + * @brief dB_{0,d}(c) = B_{0,d}(bMax, c) - B_{0,d}(0, c) + * + * @param d index of the base case + * @param u 1 + 2 * bMax^2 + * @param eULeading exp(-c / (2u)), or expm1(-c / (2u)) for the reduced form + * @param eU exp(-c / (2u)) + * @param e0 exp(-c / 2) + * @param dEi Ei(-c / (2u)) - Ei(-c / 2) + * @param c squared norm of the sample, must be > 0 + * @return dB_{0,d}(c), less u/4 for d == 0 if eULeading is the expm1 form + */ +inline double lcd_delta_b0(size_t d, double u, double eULeading, double eU, + double e0, double dEi, double c) { + assert(c > 0.00); + + if (d == 0) return 0.25 * u * eULeading - 0.25 * e0 + 0.125 * c * dEi; + if (d == 1) return -0.25 * dEi; + + double coefficient = 0.50; + double cPower = c; // c^(d-j+1), j = d + double uPower = 1.00; // u^(j-2), j = d ... built downwards below + for (size_t j = d; j > 2; --j) uPower *= u; + + double sum = 0.00; + for (size_t j = d;; --j) { + sum += coefficient * (eU / (cPower * uPower) - e0 / cPower); + if (j == 2) break; + coefficient *= static_cast(j - 2) * 2.00; + cPower *= c; + uPower /= u; + } + + return sum; +} + +/** + * @brief dBkk(k, c) = B_{k,k}(bMax, c) - B_{k,k}(0, c), the D2 term + * + * D2 = pow(2, k) * sum_i w_i * lcd_delta_bkk(k, bMax, c_i), with + * c_i = ||s_i||^2. + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param bMax upper integration bound, must be > 0 + * @param c squared norm of the sample, must be >= 0 + * @return dBkk(k, c) + * @note standard-normal target and even N only + */ +inline double lcd_delta_bkk(size_t k, double bMax, double c); + +/** + * @brief dBkk1(k, c) = B_{k,k+1}(bMax, c) - B_{k,k+1}(0, c), the D2 gradient + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param bMax upper integration bound, must be > 0 + * @param c squared norm of the sample, must be >= 0 + * @return dBkk1(k, c) + * @note standard-normal target and even N only + */ +inline double lcd_delta_bkk1(size_t k, double bMax, double c); + +/** + * @brief closed form of dBkk / dBkk1 + * + * Exposed separately so the threshold calibration test can measure exactly + * where this loses accuracy. Production callers want lcd_delta_bkk() / + * lcd_delta_bkk1(), which apply the guard. + * + * @param k half the dimension, N = 2k + * @param bMax upper integration bound + * @param c squared norm of the sample, must be > 0 + * @param shift 0 for dBkk (uses dB_{0,j}), 1 for dBkk1 (uses dB_{0,j+1}) + * @return pow(2,-k) * sum_{j=0..k} (-1)^j * binom(k,j) * dB_{0,j+shift}(c) + */ +inline double lcd_delta_bkk_closed(size_t k, double bMax, double c, + size_t shift, bool reduced = false) { + assert(k >= 1); + assert(bMax > 0.00); + assert(c > 0.00); + assert(!reduced || shift == 0); // no reduction is defined for dBkk1 + + const double u = 1.00 + 2.00 * bMax * bMax; + const double exponent = -c / (2.00 * u); + const double eU = std::exp(exponent); + const double eULeading = reduced ? std::expm1(exponent) : eU; + const double e0 = std::exp(-c / 2.00); + const double dEi = lcd_ei(exponent) - lcd_ei(-c / 2.00); + + double sum = 0.00; + double binomial = 1.00; // binom(k, j) + for (size_t j = 0; j <= k; ++j) { + const double term = + lcd_delta_b0(j + shift, u, eULeading, eU, e0, dEi, c); + sum += (j % 2 == 0) ? binomial * term : -binomial * term; + binomial *= static_cast(k - j) / static_cast(j + 1); + } + + return std::ldexp(sum, -static_cast(k)); +} + +/** + * @brief dBkk(k, 0) = integral{0 to bMax} b^(2k+1) / (1 + 2b^2)^k db + * + * The c = 0 value of the attraction term, in elementary closed form. With + * u = 1 + 2b^2 the substitution gives + * + * Z = 2^-(k+2) * integral{1 to U} (1 - 1/u)^k du + * + * = 2^-(k+2) * [ (U-1) - k*ln(U) + * + sum_{j=2..k} (-1)^j C(k,j) (U^(1-j) - 1)/(1-j) ] + * + * This is the constant lcd_delta_bkk_reduced() removes. It is O(bMax^2) by + * construction, so it belongs in the reported + * offset, never inside the optimizer's objective. + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param bMax upper integration bound, must be > 0 + * @return dBkk(k, 0) + */ +inline double lcd_delta_bkk_zero(size_t k, double bMax) { + assert(k >= 1); + assert(bMax > 0.00); + + const double u = 1.00 + 2.00 * bMax * bMax; + + double sum = (u - 1.00) - static_cast(k) * std::log(u); + + double binomial = static_cast(k); // binom(k, 1), stepped up below + for (size_t j = 2; j <= k; ++j) { + binomial *= static_cast(k - j + 1) / static_cast(j); + const double dj = static_cast(j); + const double term = + binomial * (std::pow(u, 1.00 - dj) - 1.00) / (1.00 - dj); + sum += (j % 2 == 0) ? term : -term; + } + + return std::ldexp(sum, -static_cast(k) - 2); +} + +/** + * @brief lcd_delta_bkk_zero(k, bMax) - pow(2, -k) * U / 4 + * + * The same expression with (U - 1) replaced by -1, i.e. with the O(bMax^2) + * leading term removed. The result is O(ln U), so subtracting it from the + * expm1 form of the closed dBkk introduces no cancellation. + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param bMax upper integration bound, must be > 0 + * @return dBkk(k, 0) - pow(2, -k) * (1 + 2 bMax^2) / 4 + */ +inline double lcd_delta_bkk_zero_minus_leading(size_t k, double bMax) { + assert(k >= 1); + assert(bMax > 0.00); + + const double u = 1.00 + 2.00 * bMax * bMax; + + double sum = -1.00 - static_cast(k) * std::log(u); + + double binomial = static_cast(k); + for (size_t j = 2; j <= k; ++j) { + binomial *= static_cast(k - j + 1) / static_cast(j); + const double dj = static_cast(j); + const double term = + binomial * (std::pow(u, 1.00 - dj) - 1.00) / (1.00 - dj); + sum += (j % 2 == 0) ? term : -term; + } + + return std::ldexp(sum, -static_cast(k) - 2); +} + +/** + * @brief dBkk(k, c) - dBkk(k, 0), the cancellation-free attraction term + * + * Equal to integral{0 to bMax} b^(2k+1)/(1+2b^2)^k * expm1(-c/(2(1+2b^2))) db. + * + * The constant removed per sample is lcd_delta_bkk_zero(k, bMax); D2 as a + * whole loses pow(2, k) * (sum_i w_i) * lcd_delta_bkk_zero(k, bMax). + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param bMax upper integration bound, must be > 0 + * @param c squared norm of the sample, must be >= 0 + * @return dBkk(k, c) - dBkk(k, 0) + * @note standard-normal target and even N only + */ +inline double lcd_delta_bkk_reduced(size_t k, double bMax, double c) { + assert(k >= 1); + assert(bMax > 0.00); + assert(c >= 0.00); + + if (c <= 0.00) return 0.00; // dBkk(0) - dBkk(0) + + // the expm1 form has already dropped pow(2,-k) * U/4; take off what is left + // of dBkk(k, 0), which is only O(ln U) + return lcd_delta_bkk_closed(k, bMax, c, 0, true) - + lcd_delta_bkk_zero_minus_leading(k, bMax); +} + +/** + * @brief shared body of lcd_delta_bkk / lcd_delta_bkk1 + * + * @param k half the dimension, N = 2k + * @param bMax upper integration bound + * @param c squared norm of the sample + * @param shift 0 for dBkk (uses dB_{0,j}), 1 for dBkk1 (uses dB_{0,j+1}) + * @return dBkk (shift 0) or dBkk1 (shift 1) + */ +inline double lcd_delta_bkk_impl(size_t k, double bMax, double c, + size_t shift) { + assert(k >= 1); + assert(bMax > 0.00); + assert(c >= 0.00); + + return lcd_delta_bkk_closed(k, bMax, c, shift); +} + +inline double lcd_delta_bkk(size_t k, double bMax, double c) { + return lcd_delta_bkk_impl(k, bMax, c, 0); +} + +inline double lcd_delta_bkk1(size_t k, double bMax, double c) { + return lcd_delta_bkk_impl(k, bMax, c, 1); +} + + +/** + * @brief the constant that turns the reduced objective back into the true + * distance + * + * The optimizer minimises a reduced objective in which every O(bMax^2) piece + * has been removed analytically: + * + * f_reduced = -2 * pow(2,k) * sum_i w_i * lcd_delta_bkk_reduced(k, bMax, c_i) + * + sum_ij w_i w_j * lcd_c_repulsion_reduced(bMax, T_ij) + * + * and the true distance is D = f_reduced + K, with + * + * K = A_N(bMax) + * - pow(2, k+1) * W * lcd_delta_bkk_zero(k, bMax) + * + W^2 * bMax^2 / 2, W = sum_i w_i + * + * collecting the D1 term, the per-sample attraction constant and the + * per-pair repulsion constant respectively. + * + * @param N dimension, must be even and >= 2 + * @param bMax upper integration bound, must be > 0 + * @param weightSum sum of the Dirac weights, normally 1 + * @return K + */ +inline double lcd_reported_offset(size_t N, double bMax, double weightSum) { + assert(N >= 2); + assert(N % 2 == 0); + assert(bMax > 0.00); + + const size_t k = N / 2; + const double bMaxSqrd = bMax * bMax; + + if (N == 2) { + // K = W^2 bMax^2/2 + bMax^2/2 - W bMax^2 - ln(1+bMax^2)/2 + // + W ln(1+2bMax^2)/2 + // + // Two rearrangements: + // - the three bMax^2 terms collapse to (1-W)^2 bMax^2 / 2; + // - W ln(1+2b^2)/2 - ln(1+b^2)/2 is split as + // W/2 * ln((1+2b^2)/(1+b^2)) + (W-1)/2 * ln(1+b^2), whose first factor + // is log1p(b^2/(1+b^2)). + const double weightDefect = 1.00 - weightSum; + return 0.50 * bMaxSqrd * weightDefect * weightDefect + + 0.50 * weightSum * std::log1p(bMaxSqrd / (1.00 + bMaxSqrd)) - + 0.50 * weightDefect * std::log1p(bMaxSqrd); + } + + return lcd_a_n(N, bMax) - + std::ldexp(weightSum * lcd_delta_bkk_zero(k, bMax), + static_cast(k) + 1) + + 0.50 * weightSum * weightSum * bMaxSqrd; +} + +#endif // LCD_EVEN_CLOSED_FORM_H diff --git a/plots/CMakeLists.txt b/plots/CMakeLists.txt index e15ddad..cecf92e 100644 --- a/plots/CMakeLists.txt +++ b/plots/CMakeLists.txt @@ -3,3 +3,7 @@ add_executable(generate_samples generate_samples.cpp) target_link_libraries(generate_samples PRIVATE approxLCD_static) + +add_executable(generate_samples_even generate_samples_even.cpp) + +target_link_libraries(generate_samples_even PRIVATE approxLCD_static) \ No newline at end of file diff --git a/plots/generate_samples_even.cpp b/plots/generate_samples_even.cpp new file mode 100644 index 0000000..59c96a1 --- /dev/null +++ b/plots/generate_samples_even.cpp @@ -0,0 +1,75 @@ +#include +#include +#include +#include +#include + +#include "gm_to_dirac_even_closed_form.h" + +int main(int argc, char** argv) { + if (std::getenv("GSL_RNG_SEED") == nullptr) { +#ifdef _WIN32 + _putenv_s("GSL_RNG_SEED", "42"); +#else + setenv("GSL_RNG_SEED", "42", 1); +#endif + } + + if (argc < 5 || argc > 6) { + std::cerr << "usage: " << argv[0] + << " [bMax]\n" + << " N dimension, must be even and >= 2\n" + << " L number of Dirac components\n" + << " sigma isotropic standard deviation (1 = standard normal)\n" + << " bMax integration bound, any value > 0 (default 100)\n"; + return 1; + } + + const size_t N = static_cast(std::stoi(argv[1])); + const size_t L = static_cast(std::stoi(argv[2])); + const char* outPath = argv[3]; + const double sigma = std::strtod(argv[4], nullptr); + const double bMax = (argc == 6) ? std::strtod(argv[5], nullptr) : 100.00; + + if (N == 0 || N % 2 != 0) { + std::cerr << "error: N must be even and > 0 (got " << N + << "); the closed form has no odd-N case\n"; + return 1; + } + if (L == 0 || !(sigma > 0.00) || !(bMax > 0.00)) { + std::cerr << "error: L, sigma and bMax must all be > 0\n"; + return 1; + } + + std::vector x(L * N, 0.00); + GslminimizerResult result; + gm_to_dirac_even_closed_form approx; + + // sigma == 1 reduces to the plain standard-normal solve + const bool ok = approx.approximate_isotropic(L, N, sigma, bMax, x.data(), + nullptr, &result, {}); + if (!ok) { + std::cerr << "error: LCD approximation did not converge for N=" << N + << ", L=" << L << ", sigma=" << sigma << ", bMax=" << bMax + << "\n"; + return 1; + } + + std::ofstream out(outPath); + if (!out) { + std::cerr << "error: cannot open '" << outPath << "' for writing\n"; + return 1; + } + out << std::setprecision(17); + for (size_t i = 0; i < L; ++i) { + for (size_t k = 0; k < N; ++k) { + if (k != 0) out << ','; + out << x[i * N + k]; + } + out << '\n'; + } + + std::cerr << "converged in " << result.iterations << " iterations (" + << result.elapsedTimeMicro << " us)\n"; + return 0; +} diff --git a/plots/plot_samples_even.py b/plots/plot_samples_even.py new file mode 100644 index 0000000..0ea9086 --- /dev/null +++ b/plots/plot_samples_even.py @@ -0,0 +1,253 @@ +import argparse +import os +import subprocess +import sys +import tempfile + +import numpy as np + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +# plotting constants +GRID_RES = 400 +SEED = "42" +BMAX = 5.0 + +# Figure A: standard normal approximated with more and more samples. +STANDARD_NORMAL_CASES = [ + {"sigma": 1.0, "L": 100}, + {"sigma": 1.0, "L": 200}, + {"sigma": 1.0, "L": 400}, +] + +# Figure B: isotropic sigma^2 * I at a fixed L +ISOTROPIC_CASES = [ + {"sigma": 0.5, "L": 100}, + {"sigma": 1.0, "L": 100}, + {"sigma": 2.0, "L": 100}, +] + +# Figure C: the two implementations side by side on the same problem. +COMPARISON_L = 400 + + +def resolve_generator(path): + for candidate in (path, path + ".exe"): + if os.path.exists(candidate): + return os.path.abspath(candidate) + sys.exit(f"error: generator not found at '{path}' (or '{path}.exe'). ") + + +def run_even_generator(generator, sigma, n_samples, n_dim=2): + env = dict(os.environ, GSL_RNG_SEED=SEED) + with tempfile.TemporaryDirectory() as tmp: + out_csv = os.path.join(tmp, "samples.csv") + cmd = [ + generator, + str(n_dim), + str(n_samples), + out_csv, + "%.17g" % sigma, + "%.17g" % BMAX, + ] + subprocess.run(cmd, env=env, check=True) + samples = np.loadtxt(out_csv, delimiter=",") + return samples.reshape(n_samples, n_dim) + + +def run_short_generator(generator, sigma, n_samples, n_dim=2): + """Invoke the existing quadrature sampler, for the comparison figure.""" + env = dict(os.environ, GSL_RNG_SEED=SEED) + with tempfile.TemporaryDirectory() as tmp: + out_csv = os.path.join(tmp, "samples.csv") + cmd = [generator, str(n_dim), str(n_samples), out_csv] + cmd += ["%.17g" % sigma] * n_dim + cmd += [str(int(BMAX))] + subprocess.run(cmd, env=env, check=True) + samples = np.loadtxt(out_csv, delimiter=",") + return samples.reshape(n_samples, n_dim) + + +def gaussian_density(sigma, axis_limit): + axis = np.linspace(-axis_limit, axis_limit, GRID_RES) + xx, yy = np.meshgrid(axis, axis) + quad = (xx * xx + yy * yy) / (sigma * sigma) + norm = 1.0 / (2.0 * np.pi * sigma * sigma) + return norm * np.exp(-0.5 * quad) + + +def draw_panel(ax, points, sigma, title, axis_limit=None): + if axis_limit is None: + axis_limit = 4.0 * sigma + + ax.imshow( + gaussian_density(sigma, axis_limit), + origin="lower", + extent=[-axis_limit, axis_limit, -axis_limit, axis_limit], + vmin=0.0, + aspect="equal", + ) + ax.scatter( + points[:, 0], + points[:, 1], + s=10, + c="#e8000b", + edgecolors="white", + linewidths=0.6, + zorder=3, + ) + + ax.set_xlim(-axis_limit, axis_limit) + ax.set_ylim(-axis_limit, axis_limit) + ticks = np.linspace(-axis_limit, axis_limit, 5) + ax.set_xticks(ticks) + ax.set_yticks(ticks) + ax.set_title(title, fontsize=11) + + +def build_standard_normal_figure(generator, out_path): + fig, axes = plt.subplots( + 1, len(STANDARD_NORMAL_CASES), figsize=(4.2 * len(STANDARD_NORMAL_CASES), 4.8) + ) + fig.patch.set_facecolor("white") + for ax, case in zip(np.atleast_1d(axes), STANDARD_NORMAL_CASES): + points = run_even_generator(generator, case["sigma"], case["L"]) + draw_panel(ax, points, case["sigma"], "L = %d\nΣ = I" % case["L"]) + fig.suptitle( + "Closed-form LCD (N = 2): standard normal with increasing sample counts", + fontsize=13, + y=0.99, + ) + fig.tight_layout(rect=(0, 0, 1, 0.96)) + fig.savefig(out_path, dpi=150, facecolor=fig.get_facecolor()) + plt.close(fig) + print(f"wrote {out_path}") + + +def build_isotropic_figure(generator, out_path): + fig, axes = plt.subplots( + 1, len(ISOTROPIC_CASES), figsize=(4.2 * len(ISOTROPIC_CASES), 4.8) + ) + fig.patch.set_facecolor("white") + + axis_limit = 3.0 * max(case["sigma"] for case in ISOTROPIC_CASES) + + for ax, case in zip(np.atleast_1d(axes), ISOTROPIC_CASES): + sigma = case["sigma"] + points = run_even_generator(generator, sigma, case["L"]) + radius = np.hypot(points[:, 0], points[:, 1]).max() + draw_panel( + ax, + points, + sigma, + "L = %d, σ = %g\nΣ = %g · I" % (case["L"], sigma, sigma * sigma), + axis_limit=axis_limit, + ) + fig.suptitle( + "Closed-form LCD (N = 2): isotropic Σ = σ² I", + fontsize=13, + y=0.99, + ) + fig.tight_layout(rect=(0, 0, 1, 0.96)) + fig.savefig(out_path, dpi=150, facecolor=fig.get_facecolor()) + + out_path = out_path.replace(".png", "_scaled_axes.png") + fig, axes = plt.subplots( + 1, len(ISOTROPIC_CASES), figsize=(4.2 * len(ISOTROPIC_CASES), 4.8) + ) + fig.patch.set_facecolor("white") + + for ax, case in zip(np.atleast_1d(axes), ISOTROPIC_CASES): + sigma = case["sigma"] + points = run_even_generator(generator, sigma, case["L"]) + # radius = np.hypot(points[:, 0], points[:, 1]).max() + draw_panel( + ax, + points, + sigma, + "L = %d, σ = %g\nΣ = %g · I" % (case["L"], sigma, sigma * sigma), + ) + fig.suptitle( + "Closed-form LCD (N = 2): isotropic Σ = σ² I", + fontsize=13, + y=0.99, + ) + fig.tight_layout(rect=(0, 0, 1, 0.96)) + fig.savefig(out_path, dpi=150, facecolor=fig.get_facecolor()) + + plt.close(fig) + print(f"wrote {out_path}") + + +def build_comparison_figure(even_generator, short_generator, out_path): + fig, axes = plt.subplots(1, 2, figsize=(8.4, 4.8)) + fig.patch.set_facecolor("white") + + closed_form = run_even_generator(even_generator, 1.0, COMPARISON_L) + quadrature = run_short_generator(short_generator, 1.0, COMPARISON_L) + + for ax, points, label in ( + (axes[0], closed_form, "closed form (new)"), + (axes[1], quadrature, "quadrature (existing)"), + ): + covariance = np.cov(points.T, bias=True) + draw_panel( + ax, + points, + 1.0, + "%s\ntr(Cov)/2 = %.4f" % (label, np.trace(covariance) / 2.0), + ) + + fig.suptitle( + "Closed-form vs quadrature LCD (N = 2, L = %d, bMax = %g, same seed)" + % (COMPARISON_L, BMAX), + fontsize=13, + y=0.99, + ) + fig.tight_layout(rect=(0, 0, 1, 0.96)) + fig.savefig(out_path, dpi=150, facecolor=fig.get_facecolor()) + plt.close(fig) + print(f"wrote {out_path}") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--generator", + default=os.path.join("build", "plots", "generate_samples_even"), + help="path to the compiled generate_samples_even binary", + ) + parser.add_argument( + "--short-generator", + default=os.path.join("build", "plots", "generate_samples"), + help="path to the compiled generate_samples binary (comparison figure)", + ) + parser.add_argument( + "--out-dir", + default=os.path.join("doxygen", "images"), + help="directory to write the PNG figures into", + ) + args = parser.parse_args() + + generator = resolve_generator(args.generator) + short_generator = resolve_generator(args.short_generator) + os.makedirs(args.out_dir, exist_ok=True) + + build_standard_normal_figure( + generator, os.path.join(args.out_dir, "samples_even_standard_normal.png") + ) + build_isotropic_figure( + generator, os.path.join(args.out_dir, "samples_even_isotropic.png") + ) + build_comparison_figure( + generator, + short_generator, + os.path.join(args.out_dir, "samples_even_vs_quadrature.png"), + ) + + +if __name__ == "__main__": + main() From 4d0cd0e8357b3c43333726b73288dc97f97eb3ed Mon Sep 17 00:00:00 2001 From: Balaji S B Date: Tue, 1 Sep 2026 19:25:00 +0530 Subject: [PATCH 2/2] Minor modifications --- doxygen/doxygen_gm2d_even.h | 20 ++++++------- lib/gm_to_dirac_even/lcd_even_closed_form.h | 31 ++------------------- 2 files changed, 12 insertions(+), 39 deletions(-) diff --git a/doxygen/doxygen_gm2d_even.h b/doxygen/doxygen_gm2d_even.h index 543ae46..b3986ac 100644 --- a/doxygen/doxygen_gm2d_even.h +++ b/doxygen/doxygen_gm2d_even.h @@ -63,11 +63,11 @@ * * Common parameters: * - * - L → number of Dirac components - * - N → dimension, must be even - * - bMax → integration bound, any value > 0 - * - x → initial guess and output locations (L × N) - * - wX → weights of the Dirac mixture (optional) + * - L -> number of Dirac components + * - N -> dimension, must be even + * - bMax -> integration bound, any value > 0 + * - x -> initial guess and output locations (L × N) + * - wX -> weights of the Dirac mixture (optional) * * If wX is nullptr: * @@ -136,11 +136,11 @@ * two are directly comparable. The building blocks live in * lcd_even_closed_form.h: * - * - lcd_a_n() → D1, valid for all N (odd and even) - * - lcd_delta_bkk() → the per-sample D2 term - * - lcd_delta_bkk1() → the per-sample D2 gradient term - * - lcd_c_repulsion() → the exact per-pair D3 term - * - lcd_ei() → underflow-guarded Ei + * - lcd_a_n() -> D1, valid for all N (odd and even) + * - lcd_delta_bkk() -> the per-sample D2 term + * - lcd_delta_bkk1() -> the per-sample D2 gradient term + * - lcd_c_repulsion() -> the exact per-pair D3 term + * - lcd_ei() -> Ei(x) * * \section gm_even_notes Notes * diff --git a/lib/gm_to_dirac_even/lcd_even_closed_form.h b/lib/gm_to_dirac_even/lcd_even_closed_form.h index 8d161b3..8f63956 100644 --- a/lib/gm_to_dirac_even/lcd_even_closed_form.h +++ b/lib/gm_to_dirac_even/lcd_even_closed_form.h @@ -8,8 +8,6 @@ #include #include -#include "gsl_quadrature_adaptive_gauss_kronrod.h" - /** * @file lcd_even_closed_form.h * @brief quadrature-free closed forms for the LCD / modified Cramer-von-Mises @@ -35,7 +33,7 @@ * @return Ei(x), or exactly 0.0 in case of underflow, domain error or overflow */ inline double lcd_ei(double x) { - assert(x == 0.00); + assert(x < 0.00); gsl_sf_result result; const int status = gsl_sf_expint_Ei_e(x, &result); @@ -107,27 +105,6 @@ inline double lcd_c_repulsion(double b, double c) { return 0.50 * bSqrd * std::exp(z) + 0.125 * c * lcd_ei(z); } -/** - * @brief C(b, c) - b^2/2, the repulsion term with its leading term removed - * - * Removed from D3 in total: (sum_i w_i)^2 * bMax^2 / 2. - * - * @param b upper integration bound, must be > 0 - * @param c squared distance between the two samples, must be >= 0 - * @return C(b, c) - b^2/2; exactly 0 at c = 0, so coincident samples - * (including every i == j pair) now contribute nothing - */ -inline double lcd_c_repulsion_reduced(double b, double c) { - assert(b > 0.00); - assert(c >= 0.00); - - if (c <= 0.00) return 0.00; // coincident samples - - const double bSqrd = b * b; - const double z = -c / (4.00 * bSqrd); - return 0.50 * bSqrd * std::expm1(z) + 0.125 * c * lcd_ei(z); -} - /** * @brief parameters of the raw b-integrand used by the small-c fallback */ @@ -225,10 +202,6 @@ inline double lcd_delta_bkk1(size_t k, double bMax, double c); /** * @brief closed form of dBkk / dBkk1 * - * Exposed separately so the threshold calibration test can measure exactly - * where this loses accuracy. Production callers want lcd_delta_bkk() / - * lcd_delta_bkk1(), which apply the guard. - * * @param k half the dimension, N = 2k * @param bMax upper integration bound * @param c squared norm of the sample, must be > 0 @@ -393,7 +366,7 @@ inline double lcd_delta_bkk1(size_t k, double bMax, double c) { * has been removed analytically: * * f_reduced = -2 * pow(2,k) * sum_i w_i * lcd_delta_bkk_reduced(k, bMax, c_i) - * + sum_ij w_i w_j * lcd_c_repulsion_reduced(bMax, T_ij) + * + sum_ij w_i w_j * C(bMax, T_ij) * * and the true distance is D = f_reduced + K, with *