diff --git a/.github/workflows/deploy-doxygen.yml b/.github/workflows/deploy-doxygen.yml index 894e33c..1a8d7bc 100644 --- a/.github/workflows/deploy-doxygen.yml +++ b/.github/workflows/deploy-doxygen.yml @@ -30,7 +30,7 @@ jobs: - name: Build sample generator run: | cmake -S . -B build -DCMAKE_BUILD_TYPE=Release - cmake --build build --target generate_samples -j + cmake --build build --target generate_samples generate_samples_even generate_samples_even_binf -j - name: Generate sample plots env: @@ -38,6 +38,7 @@ jobs: 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 + python plots/plot_samples_even_binf.py --generator build/plots/generate_samples_even_binf --out-dir doxygen/images - name: Build docs run: | diff --git a/doxygen/doxygen_gm2d_even_binf.h b/doxygen/doxygen_gm2d_even_binf.h new file mode 100644 index 0000000..445903b --- /dev/null +++ b/doxygen/doxygen_gm2d_even_binf.h @@ -0,0 +1,153 @@ +/** + * \page page_gm_even_binf gm_to_dirac_even_binf + * + * Gaussian-to-Dirac approximation, quadrature-free and bMax-free (even dimensions). + * + * \section gm_even_binf_overview Overview + * + * gm_to_dirac_even_binf approximates a standard normal distribution by a + * Dirac mixture with L components in N dimensions, exactly as + * gm_to_dirac_even_closed_form does, but with the limit bMax -> infinity + * taken analytically. + * + * The finite-bMax path evaluates every b-integral in closed form but still carries the integration + * bound as a parameter, and the answer it reports depends on it. Here the bound is gone: + * there is no bMax argument on any method, and ApproximateOptions::bMax is ignored. + * + * With N = 2k, \f$ c_i = \parallel s_i \parallel^2 \f$ and \f$ T_ij = \parallel s_i - s_j \parallel ^2 \f$, the limit distance + * is + * + *
\f$ D = K_k - 2 \sum_i w_i \phi_k(c_i) + \sum_{i,j, T_{ij} > 0} w_i w_j \left( \frac{T_{ij}}{8} \right) \ln\left( \frac{T_{ij}}{4} \right) \f$
+ * + * in units of \f$ \pi^{N/2} \f$, the same normalisation the other paths use, so all + * three are directly comparable. + * + * + * \section gm_even_binf_gallery Sample gallery + * + * Every point below was produced by gm_to_dirac_even_binf alone, no + * quadrature path, no finite-bMax path, no bMax anywhere. The shaded field is + * the target standard normal density and the red points are the optimized + * Dirac locations. D is the reported bMax-free mCvM distance at the optimum, + * K_k included. + * + * \image html samples_even_binf_standard_normal.png "bMax-free closed-form LCD (N = 2): standard normal at L = 100, 200, 400." width=80% + * + * Reproduce with: + * + * \code + * cmake --build build --target generate_samples_even_binf + * python plots/plot_samples_even_binf.py \ + * --generator build/plots/generate_samples_even_binf \ + * --out-dir doxygen/images + * \endcode + * + * + * \section gm_even_binf_preconditions Preconditions + * + * - N must be even and >= 2. + * - The target must be the standard normal N(0, I). There is no isotropic variant here. + * - The sample set must have zero mean, \f$ \sum_i w_i s_i = 0 \f$. + * + * + * \section gm_even_binf_zeromean The zero-mean precondition + * + * The bMax -> infinity limit exists only on the zero-mean manifold. Off it, + * the true distance diverges, + * + *
\f$ D(bMax) \sim \frac{||\mu||^2}{4} * ln(bMax^2) + const, \ \mu = \sum_i w_i s_i \f$
+ * + * and grows without bound as the integration bound is relaxed. + * + * \warning The closed form above nevertheless evaluates to a perfectly ordinary, + * finite number when it is handed a sample set with a non-zero mean. However, this result would + * not be the true distance measure, because it does not include the terms which disappeared + * when using the fact that the sample mean is zero. Thus, this result cannot be applied + * to computing the distance measure for cases without zero mean of the deterministic sample set. + * + * The class removes the hazard structurally rather than checking for it. + * The minimizer varies an unconstrained t, and every objective and gradient evaluation begins by projecting + * + *
\f$ s = t - 1 * (w^T t) \f$
+ * + * with the assembled gradient chain-ruled back by + * + *
\f$ grad_{t,q} = grad_{s,q} - w_q * \sum_i grad_{s,i} \f$
+ * + * so no iterate can leave the manifold and the finite non-distance is unreachable. + * The same projection runs inside the standalone distance and derivative entry points, + * which therefore report the value for the mean-corrected input rather than for the raw input. + * + * correctMean() after minimize() is not a substitute for this and must not be relied on as one. + * It is called at the end because it is exactly the same projection and so turns the returned + * optimizer variable t into the sample set s that the reported objective belongs to. + * + * \section gm_even_binf_interface Interface + * + * Template parameter: + * + * - T in {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(...) + * + * for raw pointers (T*), GSL vectors (gsl_vector / gsl_vector_float) and GSL + * matrices (gsl_matrix / gsl_matrix_float). x is L x N, row major. + * + * If wX is nullptr, uniform weights are assumed. + * + * \note modified_van_mises_distance_sq_derivative() returns \f$ grad_s \f$, the + * derivative with respect to the sample locations. That is not the gradient + * the minimizer works with; approximate() uses \f$ grad_t \f$, which additionally + * carries the projection chain rule above. + * + * + * \section example_gm_even_binf_raw Example (Raw Pointer) + * + * \code + * gm_to_dirac_even_binf approx; + * + * bool ok = approx.approximate( + * L, + * N, // must be even + * x, // initial guess / output (L x N), returned zero-mean + * wX, // weights (optional) + * &result, + * options // options.bMax is ignored + * ); + * \endcode + * + * + * \section example_gm_even_binf_distance Example (true distance) + * + * \code + * gm_to_dirac_even_binf approx; + * + * double distance = 0.0; + * approx.modified_van_mises_distance_sq(&distance, L, N, x, wX); + * // includes K_k, so this is the true bMax-free mCvM distance of the + * // mean-corrected sample set + * \endcode + * + * + * \section gm_even_binf_math Closed forms + * + * The building blocks live in lcd_even_binf_closed_form.h: + * + * - lcd_binf_bkk_zero() -> \f$ B_{k,k}(0, c) \f$, the attraction term + * - lcd_binf_bkk_zero_deriv() -> \f$ B_{k,k}'(0, c) \f$, its derivative + * - lcd_binf_phi() -> \f$ \phi_k(c) \f$, the per-sample potential + * - lcd_binf_phi_deriv() -> \f$ \phi_k'(c) \f$, its derivative + * - lcd_binf_constant() -> \f$ K_k \f$ + * + * Ei is taken from lcd_ei() of lcd_even_closed_form.h. + * + * K_k is a pure constant. It does not move the optimum, and is carried so + * that the reported value is the true distance; it plays the role + * lcd_reported_offset() plays on the finite-bMax path. + * + */ diff --git a/doxygen/doxygen_main_page.h b/doxygen/doxygen_main_page.h index ce1a506..28d19bb 100644 --- a/doxygen/doxygen_main_page.h +++ b/doxygen/doxygen_main_page.h @@ -42,14 +42,12 @@ * Implementations: * - \ref page_gm_short "gm_to_dirac_short" * - \ref page_gm_stddev "gm_to_dirac_short_standard_normal_deviation" - * - \ref page_gm_even "gm_to_dirac_even_closed_form" — quadrature-free; + * - \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. + * - \ref page_gm_even_binf "gm_to_dirac_even_binf" : quadrature-free AND + * bMax-free, the bMax -> infinity limit of the entry above; requires even + * N, a standard normal target and a zero-mean sample set, which it enforces by optimising on the zero-mean manifold * * \section samples Sample Gallery * diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index bbe9ec3..ead3b92 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -40,6 +40,10 @@ file(GLOB GM_TO_DIRAC_C_WRAPPER_SOURCES ${GM_TO_DIRAC_C_WRAPPER_DIR}/*.cpp) 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) +# gm_to_dirac_even_binf sources (bMax-free limit of the even-N closed form) +set(GM_TO_DIRAC_EVEN_BINF_DIR ${CMAKE_CURRENT_SOURCE_DIR}/gm_to_dirac_even_binf) +file(GLOB GM_TO_DIRAC_EVEN_BINF_SOURCES ${GM_TO_DIRAC_EVEN_BINF_DIR}/*.cpp) + # common sources set(COMMON_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/debug/capture_time.cpp @@ -69,6 +73,7 @@ add_library(approxLCD_static STATIC ${GM_TO_DIRAC_SOURCES} ${GM_TO_DIRAC_C_WRAPPER_SOURCES} ${GM_TO_DIRAC_EVEN_SOURCES} + ${GM_TO_DIRAC_EVEN_BINF_SOURCES} ) target_include_directories(approxLCD_static PUBLIC ${COMMON_INCLUDES} @@ -77,6 +82,7 @@ target_include_directories(approxLCD_static PUBLIC ${GM_TO_DIRAC_DIR} ${GM_TO_DIRAC_C_WRAPPER_DIR} ${GM_TO_DIRAC_EVEN_DIR} + ${GM_TO_DIRAC_EVEN_BINF_DIR} ) set_target_properties(approxLCD_static PROPERTIES OUTPUT_NAME "approxLCD") target_link_libraries(approxLCD_static PUBLIC GSL::gsl GSL::gslcblas OpenMP::OpenMP_CXX) @@ -89,6 +95,7 @@ add_library(approxLCD_shared SHARED ${GM_TO_DIRAC_SOURCES} ${GM_TO_DIRAC_C_WRAPPER_SOURCES} ${GM_TO_DIRAC_EVEN_SOURCES} + ${GM_TO_DIRAC_EVEN_BINF_SOURCES} ) target_include_directories(approxLCD_shared PUBLIC ${COMMON_INCLUDES} @@ -97,6 +104,7 @@ target_include_directories(approxLCD_shared PUBLIC ${GM_TO_DIRAC_DIR} ${GM_TO_DIRAC_C_WRAPPER_DIR} ${GM_TO_DIRAC_EVEN_DIR} + ${GM_TO_DIRAC_EVEN_BINF_DIR} ) # platform-specific options diff --git a/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.cpp b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.cpp new file mode 100644 index 0000000..992db76 --- /dev/null +++ b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.cpp @@ -0,0 +1,303 @@ +#include "gm_to_dirac_even_binf.h" + +#include +#include + +#include +#include +#include +#include + +#include "gsl_minimizer.h" +#include "gsl_utils_view_helper.h" +#include "gsl_utils_weight_helper.h" + +namespace { + constexpr double kBinfFailureValue = std::numeric_limits::quiet_NaN(); +} + +/** + * @brief objective and gradient of the bMax-free distance + * + * The order matters: project t onto the zero-mean manifold and refresh the + * c_i / T_ij caches first, then accumulate the gradient with respect to s. + * With chainRuleToT set, the assembled gradient is finally mapped back to the + * minimizer's variable t; without it, grad_s is returned as is. + * + * A sample projected exactly onto the origin is outside the closed form's + * domain (see GMToDiracEvenBinfOptimizationParams::hasDegenerateSample); it + * yields NaN here rather than the abort lcd_ei(0) would trigger. + * + * @param t optimizer variable (or, for the public entry points, the sample locations), L * N + * @param params optimization parameters + * @param f objective accumulator, may be nullptr + * @param grad gradient accumulator, may be nullptr + * @param chainRuleToT true for grad_t, false for grad_s + */ +template +void gm_to_dirac_even_binf::evaluate( + const gsl_vector* t, GMToDiracEvenBinfOptimizationParams* params, double* f, + gsl_vector* grad, bool chainRuleToT) { + if (f) *f = 0.00; + if (grad) gsl_vector_set_zero(grad); + + params->update(t); // s = t - 1 * (w^T t), then the caches + + if (params->hasDegenerateSample()) { + if (f) *f = kBinfFailureValue; + if (grad) + for (size_t i = 0; i < grad->size; ++i) grad->data[i] = kBinfFailureValue; + return; + } + + calculateAttraction(params, f, grad); + calculateRepulsion(params, f, grad); + + if (grad && chainRuleToT) projectGradient(params, grad); +} + +template +double gm_to_dirac_even_binf::modified_van_mises_distance_sq( + const gsl_vector* t, void* params) { + double d = 0.00; + evaluate(t, static_cast(params), &d, + nullptr, true); + return d; +} + +template +void gm_to_dirac_even_binf::modified_van_mises_distance_sq_derivative( + const gsl_vector* t, void* params, gsl_vector* grad) { + evaluate(t, static_cast(params), + nullptr, grad, true); +} + +template +void gm_to_dirac_even_binf::combined_distance_metric(const gsl_vector* t, + void* params, double* f, + gsl_vector* grad) { + evaluate(t, static_cast(params), f, + grad, true); +} + +// raw pointer + +template +bool gm_to_dirac_even_binf::approximate(size_t L, size_t N, 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, vectorViewX.get(), vectorViewWX.get(), result, + options); +} + +template +void gm_to_dirac_even_binf::modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, T* x, + const T* wX) { + GSLVectorView vectorViewX(x, L * N); + GSLVectorView vectorViewWX(wX, L); + modified_van_mises_distance_sq(distance, L, N, vectorViewX.get(), + vectorViewWX.get()); +} + +template +void gm_to_dirac_even_binf::modified_van_mises_distance_sq_derivative( + T* gradient, size_t L, size_t N, 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, + vectorViewX.get(), + vectorViewWX.get()); +} + +// gsl matrix + +template +bool gm_to_dirac_even_binf::approximate(size_t L, size_t N, 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, vectorViewX.get(), wX, result, options); +} + +template +void gm_to_dirac_even_binf::modified_van_mises_distance_sq( + T* distance, size_t L, size_t N, GSLMatrixType* x, + const GSLVectorType* wX) { + assert(x->size1 == L); + assert(x->size2 == N); + GSLVectorView vectorViewX(x); + modified_van_mises_distance_sq(distance, L, N, vectorViewX.get(), wX); +} + +template +void gm_to_dirac_even_binf::modified_van_mises_distance_sq_derivative( + GSLMatrixType* gradient, size_t L, size_t N, 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, + vectorViewX.get(), wX); +} + +// float + +template <> +bool gm_to_dirac_even_binf::approximate( + size_t L, size_t N, gsl_vector_float* x, const gsl_vector_float* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + assert(x != nullptr); + 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_binf doubleApprox; + const bool success = + doubleApprox.approximate(L, N, 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_binf::modified_van_mises_distance_sq( + float* distance, size_t L, size_t N, 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_binf doubleApprox; + doubleApprox.modified_van_mises_distance_sq( + &distanceDouble, L, N, vectorViewX.get(), vectorViewWX.get()); + *distance = static_cast(distanceDouble); +} + +template <> +void gm_to_dirac_even_binf::modified_van_mises_distance_sq_derivative( + gsl_vector_float* gradient, size_t L, size_t N, 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_binf doubleApprox; + doubleApprox.modified_van_mises_distance_sq_derivative( + gradientDouble, L, N, 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_binf::approximate( + size_t L, size_t N, gsl_vector* x, const gsl_vector* wX, + GslminimizerResult* result, const ApproximateOptions& options) { + assert(x != nullptr); + assert(x->size == L * N); + if (!preconditionsHold(L, N)) 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); + GMToDiracEvenBinfOptimizationParams params(wXHelper.get(), N, L); + + // refuse a start that projects a sample onto the origin before the + // minimizer ever sees it; L == 1 always does, having only the point 0 on + // the zero-mean manifold + params.update(x); + if (params.hasDegenerateSample()) return false; + + 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); + + // x holds the optimizer variable t; this is the very projection every + // evaluation applied internally, so it returns the sample set s that the + // reported objective belongs to + correctMean(x, params.wX, L, N); + + return status == GSL_SUCCESS; +} + +template <> +void gm_to_dirac_even_binf::modified_van_mises_distance_sq( + double* distance, size_t L, size_t N, gsl_vector* x, const gsl_vector* wX) { + assert(distance != nullptr); + if (!preconditionsHold(L, N)) { + *distance = kBinfFailureValue; + return; + } + + GSLWeightHelper wXHelper(wX, L); + GMToDiracEvenBinfOptimizationParams optiParams(wXHelper.get(), N, L); + + evaluate(x, &optiParams, distance, nullptr, false); +} + +template <> +void gm_to_dirac_even_binf::modified_van_mises_distance_sq_derivative( + gsl_vector* gradient, size_t L, size_t N, gsl_vector* x, + const gsl_vector* wX) { + assert(gradient != nullptr); + if (!preconditionsHold(L, N)) { + for (size_t i = 0; i < gradient->size; ++i) + gradient->data[i] = kBinfFailureValue; + return; + } + + GSLWeightHelper wXHelper(wX, L); + GMToDiracEvenBinfOptimizationParams optiParams(wXHelper.get(), N, L); + + // grad_s, the gradient w.r.t. the sample locations; the minimizer's grad_t + // additionally carries the projection chain rule + evaluate(x, &optiParams, nullptr, gradient, false); +} + +template class gm_to_dirac_even_binf; +template class gm_to_dirac_even_binf; diff --git a/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.h b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.h new file mode 100644 index 0000000..3a01d88 --- /dev/null +++ b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.h @@ -0,0 +1,287 @@ +#ifndef GM_TO_DIRAC_EVEN_BINF_H +#define GM_TO_DIRAC_EVEN_BINF_H + +#include "approximate_options.h" +#include "gm_to_dirac_even_binf_optimization_params.h" +#include "gsl_minimizer.h" +#include "gsl_vector_matrix_types.h" + +/** + * @brief bMax-free Gaussian-to-Dirac approximation for an even-dimensional standard normal target + * + * Same objective family as gm_to_dirac_even_closed_form, but the limit bMax -> infinity has been + * taken analytically (see lcd_even_binf_closed_form.h), so no integration bound appears anywhere in the API. + * + * @note PRECONDITIONS + * + * - N must be EVEN and >= 2 and L >= 1. approximate() returns false and the void-returning entry points write NaN + * to their output rather than leaving a stale buffer. + * - The target must be the STANDARD normal N(0, I). + * - The sample set must have ZERO MEAN, sum_i w_i s_i = 0. The limit exists only there; for a non-zero mean the + * closed form still evaluates to a finite number that is NOT the distance. + * + * The minimizer optimises over an unconstrained t and every objective and + * gradient evaluation begins by projecting + * + * s = t - 1 * (w^T t) + * + * with the assembled gradient chain-ruled back by + * + * grad_{t,q} = grad_{s,q} - w_q * sum_i grad_{s,i} + * + * so no iterate can ever leave the zero-mean manifold and the finite non-distance is unreachable. + * The same projection runs inside the standalone distance and derivative entry points, which therefore report + * the value for the MEAN-CORRECTED input rather than for the raw input. + * + * @note DEGENERATE SAMPLE SETS. phi_k is undefined at c = 0 and evaluating it there reaches lcd_ei(0), + * which GSL's default error handler turns into a process abort. Any input whose projection puts a sample + * exactly on the origin is therefore refused the same way odd N is: approximate() returns false before + * the minimizer starts and the void entry points write NaN. + * That covers L == 1, whose only zero-mean sample set is the single point 0 and inputs such as L identical locations. + * + * @tparam T float or double; computation is always performed in double and cast at the interface boundary + */ +template +class gm_to_dirac_even_binf { + public: + using GSLVectorType = typename GSLTemplateTypeAlias::VectorType; + using GSLMatrixType = typename GSLTemplateTypeAlias::MatrixType; + + // clang-format off + /** + * @brief approximate a standard normal by L Dirac components + * + * On return x holds the projected, zero-mean sample set s, not the raw + * optimizer variable t. + * + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @param x initial guess and output locations (L * N) + * @param wX weights, nullptr for uniform + * @param result minimizer result + * @param options minimizer options; ApproximateOptions::bMax is IGNORED, + * there is no integration bound on this path + * @return true on success, false otherwise (including odd N) + */ + bool approximate(size_t L, + size_t N, + T* x, + const T* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}); + + /** + * @brief true bMax-free mCvM distance, K_k included + * + * Reports the distance of the mean-corrected sample set + * + * @param distance output distance; set to NaN if the preconditions fail + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @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, + T* x, + const T* wX = nullptr); + + /** + * @brief gradient of the true bMax-free mCvM distance w.r.t. the SAMPLES + * + * Returns grad_s, the derivative with respect to the sample locations themselves: + * + * grad_s,q = -4 w_q phi_k'(c_q) s_q + * + 0.5 w_q sum_{j != q, T_qj > 0} w_j (s_q - s_j) + * * (ln(T_qj / 4) + 1) + * + * This is NOT the gradient the minimizer works with. approximate() varies + * the unconstrained t behind s = t - 1 * (w^T t) and therefore uses + * grad_t,q = grad_s,q - w_q * sum_i grad_s,i. The two agree only up to that + * projection; grad_s is exposed here because it is the derivative of the + * distance as a function of the sample set, which is what a caller asking + * for "the gradient at x" means. + * + * @param gradient output gradient (L * N); filled with NaN if the + * preconditions fail + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @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, + T* x, + const T* wX = nullptr); + + /** + * @brief approximate a standard normal by L Dirac components + * + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @param x initial guess and output locations (L * N) + * @param wX weights, nullptr for uniform + * @param result minimizer result + * @param options minimizer options; ApproximateOptions::bMax is IGNORED + * @return true on success, false otherwise (including odd N) + */ + bool approximate(size_t L, + size_t N, + GSLVectorType* x, + const GSLVectorType* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}); + + /** + * @brief true bMax-free mCvM distance, K_k included + * + * @param distance output distance; set to NaN if the preconditions fail + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @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, + GSLVectorType* x, + const GSLVectorType* wX = nullptr); + + /** + * @brief gradient of the true bMax-free mCvM distance w.r.t. the SAMPLES + * + * Returns grad_s, not the minimizer's grad_t; see the raw-pointer overload. + * + * @param gradient output gradient (L * N); filled with NaN if the + * preconditions fail + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @param x sample locations (L * N) + * @param wX weights, nullptr for uniform + */ + void modified_van_mises_distance_sq_derivative(GSLVectorType* gradient, + size_t L, + size_t N, + GSLVectorType* x, + const GSLVectorType* wX = nullptr); + + /** + * @brief approximate a standard normal by L Dirac components + * + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @param x initial guess and output locations (L rows, N columns) + * @param wX weights, nullptr for uniform + * @param result minimizer result + * @param options minimizer options; ApproximateOptions::bMax is IGNORED + * @return true on success, false otherwise (including odd N) + */ + bool approximate(size_t L, + size_t N, + GSLMatrixType* x, + const GSLVectorType* wX = nullptr, + GslminimizerResult* result = nullptr, + const ApproximateOptions& options = ApproximateOptions{}); + + /** + * @brief true bMax-free mCvM distance, K_k included + * + * @param distance output distance; set to NaN if the preconditions fail + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @param x sample locations (L rows, N columns) + * @param wX weights, nullptr for uniform + */ + void modified_van_mises_distance_sq(T* distance, + size_t L, + size_t N, + GSLMatrixType* x, + const GSLVectorType* wX = nullptr); + + /** + * @brief gradient of the true bMax-free mCvM distance w.r.t. the SAMPLES + * + * Returns grad_s, not the minimizer's grad_t; see the raw-pointer overload. + * + * @param gradient output gradient (L rows, N columns); filled with NaN if + * the preconditions fail + * @param L number of Dirac components, must be >= 1 + * @param N dimension, must be even and >= 2 + * @param x sample locations (L rows, N columns) + * @param wX weights, nullptr for uniform + */ + void modified_van_mises_distance_sq_derivative(GSLMatrixType* gradient, + size_t L, + size_t N, + GSLMatrixType* x, + const GSLVectorType* wX = nullptr); + // clang-format on + + private: + /// shared evaluation body; chainRuleToT selects grad_t (the minimizer's + /// variable) over grad_s (the sample locations) + static void evaluate(const gsl_vector* t, + GMToDiracEvenBinfOptimizationParams* params, double* f, + gsl_vector* grad, bool chainRuleToT); + + // gsl_minimizer entry points; these always work in t + static double modified_van_mises_distance_sq(const gsl_vector* t, + void* params); + static void modified_van_mises_distance_sq_derivative(const gsl_vector* t, + void* params, + gsl_vector* grad); + static void combined_distance_metric(const gsl_vector* t, void* params, + double* f, gsl_vector* grad); + + static inline void calculateAttraction( + GMToDiracEvenBinfOptimizationParams* params, double* f, gsl_vector* grad); + + static inline void calculateRepulsion( + GMToDiracEvenBinfOptimizationParams* params, double* f, gsl_vector* grad); + + static inline void projectGradient( + const GMToDiracEvenBinfOptimizationParams* params, 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); +}; + +#include "gm_to_dirac_even_binf.tpp" + +template <> +bool gm_to_dirac_even_binf::approximate( + size_t L, size_t N, gsl_vector_float* x, const gsl_vector_float* wX, + GslminimizerResult* result, const ApproximateOptions& options); + +template <> +bool gm_to_dirac_even_binf::approximate( + size_t L, size_t N, gsl_vector* x, const gsl_vector* wX, + GslminimizerResult* result, const ApproximateOptions& options); + +template <> +void gm_to_dirac_even_binf::modified_van_mises_distance_sq( + float* distance, size_t L, size_t N, gsl_vector_float* x, + const gsl_vector_float* wX); + +template <> +void gm_to_dirac_even_binf::modified_van_mises_distance_sq( + double* distance, size_t L, size_t N, gsl_vector* x, const gsl_vector* wX); + +template <> +void gm_to_dirac_even_binf::modified_van_mises_distance_sq_derivative( + gsl_vector_float* gradient, size_t L, size_t N, gsl_vector_float* x, + const gsl_vector_float* wX); + +template <> +void gm_to_dirac_even_binf::modified_van_mises_distance_sq_derivative( + gsl_vector* gradient, size_t L, size_t N, gsl_vector* x, + const gsl_vector* wX); + +extern template class gm_to_dirac_even_binf; +extern template class gm_to_dirac_even_binf; + +#endif // GM_TO_DIRAC_EVEN_BINF_H diff --git a/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.tpp b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.tpp new file mode 100644 index 0000000..3d523b4 --- /dev/null +++ b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf.tpp @@ -0,0 +1,175 @@ +#ifndef GM_TO_DIRAC_EVEN_BINF_TPP +#define GM_TO_DIRAC_EVEN_BINF_TPP + +#include +#include +#include + +/** + * @brief check the preconditions of the bMax-free even-N closed-form path + * + * There is no bMax to validate; what remains is the parity of N and the + * component count. + * + * @param L number of Dirac components + * @param N dimension + * @return true if the closed form applies + */ +template +inline bool gm_to_dirac_even_binf::preconditionsHold(size_t L, size_t N) { + assert(N % 2 == 0); + assert(N >= 2); + assert(L >= 1); + + return N >= 2 && N % 2 == 0 && L >= 1; +} + +/** + * @brief attraction term: contributes K_k - 2 * sum_i w_i phi_k(c_i) to f and + * its gradient to grad + * + * f += K_k - 2 * sum_i w_i * phi_k(c_i) + * + * grad += -4 * w_q * phi_k'(c_q) * s_q + * + * The factor 4 rather than 2 comes from dc_q / ds_q = 2 s_q. + * + * @param params optimization parameters + * @param f objective accumulator, may be nullptr + * @param grad gradient accumulator with respect to s, may be nullptr + */ +template +inline void gm_to_dirac_even_binf::calculateAttraction( + GMToDiracEvenBinfOptimizationParams* 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 gsl_vector* wX = params->wX; + + if (f) { + double sum = 0.00; + for (size_t i = 0; i < L; ++i) + sum += wX->data[i] * lcd_binf_phi(k, params->cSqrdNorm[i]); + + *f += params->constant - 2.00 * sum; + } + + if (grad) { + for (size_t q = 0; q < L; ++q) { + const double factor = + -4.00 * wX->data[q] * lcd_binf_phi_deriv(k, params->cSqrdNorm[q]); + for (size_t d = 0; d < N; ++d) + grad->data[q * N + d] += factor * params->sample(q, d); + } + } +} + +/** + * @brief repulsion term: contributes the T ln T double sum to f and its + * gradient to grad + * + * f += sum_{i,j, T_ij > 0} w_i * w_j * (T_ij / 8) * ln(T_ij / 4) + * + * grad += 0.5 * w_q * sum_{j != q, T_qj > 0} w_j * (s_q - s_j) + * * (ln(T_qj / 4) + 1) + * + * The 1/2 rather than 1/4 is because the double sum reaches index q both as + * i and as j. Coincident samples (T <= 0) are skipped, exactly as + * gm_to_dirac_even_closed_form::calculateD3 does; the term and its gradient + * both tend to zero there. + * + * @param params optimization parameters, caches must be up to date + * @param f objective accumulator, may be nullptr + * @param grad gradient accumulator with respect to s, may be nullptr + */ +template +inline void gm_to_dirac_even_binf::calculateRepulsion( + GMToDiracEvenBinfOptimizationParams* params, double* f, gsl_vector* grad) { + const size_t L = params->L; + const size_t N = params->N; + const gsl_vector* wX = params->wX; + + double repulsion = 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 localDistSq = params->distanceSq(i, j); + + // coincident samples + if (localDistSq <= 0.00) continue; + + const double wXiwXj = wXi * wX->data[j]; + const double logTerm = std::log(0.25 * localDistSq); + + if (f) repulsion += wXiwXj * 0.125 * localDistSq * logTerm; + + if (!grad) continue; + + const double constFactor = 0.50 * wXiwXj * (logTerm + 1.00); + for (size_t d = 0; d < N; ++d) + grad->data[i * N + d] += + constFactor * (params->sample(i, d) - params->sample(j, d)); + } + } + + if (f) *f += repulsion; +} + +/** + * @brief chain-rule the gradient back through the zero-mean projection + * + * With s = t - 1 * (w^T t), + * + * grad_{t,q} = grad_{s,q} - w_q * sum_i grad_{s,i} + * + * Applied in place to the assembled ds-gradient. + * + * @param params optimization parameters + * @param grad gradient with respect to s on entry, with respect to t on exit + */ +template +inline void gm_to_dirac_even_binf::projectGradient( + const GMToDiracEvenBinfOptimizationParams* params, gsl_vector* grad) { + const size_t L = params->L; + const size_t N = params->N; + const gsl_vector* wX = params->wX; + + std::vector total(N, 0.00); + for (size_t i = 0; i < L; ++i) + for (size_t d = 0; d < N; ++d) total[d] += grad->data[i * N + d]; + + for (size_t q = 0; q < L; ++q) { + const double wXq = wX->data[q]; + for (size_t d = 0; d < N; ++d) grad->data[q * N + d] -= wXq * total[d]; + } +} + +/** + * @brief subtract the weighted mean from every sample, in place + * + * This is exactly the projection s = t - 1 * (w^T t) that every evaluation + * applies internally, so calling it once after minimize() turns the returned + * optimizer variable t into the sample set s the objective was actually evaluated at. + * + * @param x sample locations (L * N), modified in place + * @param wX weights of the Dirac mixture + * @param L number of Dirac components + * @param N dimension + */ +template +inline void gm_to_dirac_even_binf::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_BINF_TPP diff --git a/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf_optimization_params.h b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf_optimization_params.h new file mode 100644 index 0000000..1ec9f14 --- /dev/null +++ b/lib/gm_to_dirac_even_binf/gm_to_dirac_even_binf_optimization_params.h @@ -0,0 +1,152 @@ +#ifndef GM_TO_DIRAC_EVEN_BINF_OPTIMIZATION_PARAMS_H +#define GM_TO_DIRAC_EVEN_BINF_OPTIMIZATION_PARAMS_H + +#include +#include + +#include +#include + +#include "gsl_minimizer.h" +#include "lcd_even_binf_closed_form.h" +#include "squared_euclidean_distance_utils.h" + +/** + * @brief optimization parameters for the bMax-free closed-form even-N Gaussian-to-Dirac approximation + * + * Mirrors GMToDiracEvenOptimizationParams minus bMax and everything derived from it: + * there is no reportedOffset and no objectiveOffset, because the bMax-free objective carries + * its only constant, lcd_binf_constant(k), directly. + * + * It holds one thing the finite-bMax struct does not: the PROJECTED sample buffer. + * The minimizer's variable is t, but the distance is a function of + * + * s_i = t_i - sum_j w_j t_j + * + * so c_i and T_ij must be computed from s, never from t. update() performs + * that projection and then refreshes both caches from the result. + */ +struct GMToDiracEvenBinfOptimizationParams + : public GslMinimizerOptimizationParams { + public: + /** + * @brief Construct a new GMToDiracEvenBinfOptimizationParams 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 + */ + GMToDiracEvenBinfOptimizationParams(const gsl_vector* wX, size_t N, size_t L) + : GslMinimizerOptimizationParams(L, N), + wX(wX), + k(N / 2), + constant(lcd_binf_constant(N / 2)), + sProjected(L * N, 0.00), + cSqrdNorm(L, 0.00) { + assert(wX != nullptr); + assert(wX->size == L); + assert(N % 2 == 0); + assert(N >= 2); + assert(L >= 1); + + if (L > 1) + squaredEuclideanDistanceUtilLL = + new SquaredEuclideanDistance_LL_vectorized(L, N); + } + + /** + * @brief Destroy the GMToDiracEvenBinfOptimizationParams object + */ + ~GMToDiracEvenBinfOptimizationParams() { + if (squaredEuclideanDistanceUtilLL) delete squaredEuclideanDistanceUtilLL; + } + + /** + * @brief project t onto the zero-mean manifold and refresh the caches + * + * Writes s_i = t_i - sum_j w_j t_j into sProjected, then fills cSqrdNorm + * and the pairwise distances from s. Must be called whenever t changes, + * before reading sample(), cSqrdNorm or distanceSq(). + * + * @param t optimizer variable, L * N + */ + inline void update(const gsl_vector* t) { + assert(t->size == L * N); + + // s = t - 1 * (w^T t) + 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 * t->data[i * N + d]; + } + + for (size_t i = 0; i < L; ++i) { + double sum = 0.00; + for (size_t d = 0; d < N; ++d) { + const double sid = t->data[i * N + d] - mean[d]; + sProjected[i * N + d] = sid; + sum += sid * sid; + } + cSqrdNorm[i] = sum; + } + + if (!squaredEuclideanDistanceUtilLL) return; + + const gsl_matrix sMatrix = + gsl_matrix_view_array(sProjected.data(), L, N).matrix; + squaredEuclideanDistanceUtilLL->calculateDistance(&sMatrix, nullptr); + } + + /** + * @brief true if any projected sample sits exactly on the origin + * + * c_i == 0 is outside the domain of lcd_binf_phi(): it reaches lcd_ei(0), + * which aborts through GSL's default error handler. + * + * Reads the cache, so update() must have run. + * + * @return true if the closed form cannot be evaluated for this sample set + */ + inline bool hasDegenerateSample() const { + for (size_t i = 0; i < L; ++i) + if (!(cSqrdNorm[i] > 0.00)) return true; + return false; + } + + /** + * @brief cached projected sample coordinate s_{i,d} + * + * @param i index of the sample + * @param d index of the dimension + * @return s_{i,d} + */ + inline double sample(size_t i, size_t d) const { + return sProjected[i * N + d]; + } + + /** + * @brief cached squared distance T_ij between projected samples i and j + * + * The projection is a common translation of every sample, so T_ij is the + * same for t and for s; it is computed from s regardless, to keep one + * source of truth. + * + * @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 size_t k; ///< N / 2 + const double constant; ///< K_k, see lcd_binf_constant() + + std::vector sProjected; ///< s = t - 1 * (w^T t), size L * N + std::vector cSqrdNorm; ///< c_i = ||s_i||^2, size L + SquaredEuclideanDistanceUtilsLL* squaredEuclideanDistanceUtilLL = nullptr; +}; + +#endif // GM_TO_DIRAC_EVEN_BINF_OPTIMIZATION_PARAMS_H diff --git a/lib/gm_to_dirac_even_binf/lcd_even_binf_closed_form.h b/lib/gm_to_dirac_even_binf/lcd_even_binf_closed_form.h new file mode 100644 index 0000000..1a4936b --- /dev/null +++ b/lib/gm_to_dirac_even_binf/lcd_even_binf_closed_form.h @@ -0,0 +1,307 @@ +#ifndef LCD_EVEN_BINF_CLOSED_FORM_H +#define LCD_EVEN_BINF_CLOSED_FORM_H + +#include +#include +#include + +#include "lcd_even_closed_form.h" // lcd_ei() + +/** + * @file lcd_even_binf_closed_form.h + * @brief bMax-free closed forms for the LCD / modified Cramer-von-Mises distance of + * an even-dimensional standard normal target + * + * lcd_even_closed_form.h evaluates the mCvM distance for a finite integration bound bMax. + * This header takes the limit bMax -> infinity analytically, so bMax disappears from the problem entirely. + * + * With N = 2k, samples s_i, weights w_i summing to 1, c_i = ||s_i||^2 and + * T_ij = ||s_i - s_j||^2 the limit distance is + * + * D = K_k - 2 * sum_i w_i * phi_k(c_i) + * + sum_{i,j, T_ij > 0} w_i w_j * (T_ij / 8) * ln(T_ij / 4) + * + * @note PRECONDITION: the limit exists only for a ZERO-MEAN sample set, sum_i w_i s_i = 0. + * For a non-zero mean the true distance diverges like ||mu||^2 / 4 * ln(bMax^2); + * the formulas here still return a finite number in that case, but that number is NOT the distance. + * gm_to_dirac_even_binf removes the hazard by optimising on the zero-mean manifold and projecting + * before every evaluation. + * + * @note PRECONDITION: the target must be the STANDARD normal N(0, I) and + * N = 2k must be EVEN. + * + * @warning c == 0 IS A DOMAIN ERROR, not merely inaccurate: it reaches lcd_ei(0), + * whose GSL_EDOM result GSL's default error handler turns into an + * abort long before lcd_ei() can map it to 0.0. Every entry point here + * asserts c > 0, and gm_to_dirac_even_binf refuses a sample set that projects + * any component onto the origin. + */ + +/// Euler-Mascheroni constant +constexpr double lcdBinfEulerGamma = 0.5772156649015328606; + +/** + * @brief B_{0,d}(0, c), the lower-limit base case of the attraction term + * + * With e0 = exp(-c/2) and Ei0 = Ei(-c/2), + * + * B_{0,0}(0,c) = e0/4 + (c/8) * Ei0 + * + * B_{0,1}(0,c) = -Ei0/4 + * + * B_{0,d}(0,c) = e0 * sum_{m=2..d} (d-2)! * 2^(d-m-1) + * / ( (m-2)! * c^(d-m+1) ) for d >= 2 + * + * The d >= 2 sum is accumulated downwards from m = d, where the term is + * 1 / (2c), using term_{m-1} = term_m * 2 (m - 2) / c. That is the same + * recurrence lcd_delta_b0() uses and it keeps both (d-2)! and c^(d-1) out of + * the intermediate values. + * + * @param d index of the base case + * @param e0 exp(-c / 2) + * @param ei0 Ei(-c / 2) + * @param c squared norm of the sample, must be > 0 + * @return B_{0,d}(0, c) + */ +inline double lcd_binf_b0_zero(size_t d, double e0, double ei0, double c) { + assert(c > 0.00); + + if (d == 0) return 0.25 * e0 + 0.125 * c * ei0; + if (d == 1) return -0.25 * ei0; + + double term = 0.50 / c; // m = d + double sum = term; + for (size_t m = d; m > 2; --m) { + term *= 2.00 * static_cast(m - 2) / c; // step down to m - 1 + sum += term; + } + + return e0 * sum; +} + +/** + * @brief d/dc B_{0,d}(0, c) + * + * Using d/dc Ei(-c/2) = e0 / c, + * + * dB_{0,0}/dc = Ei0 / 8 + * + * dB_{0,1}/dc = -e0 / (4c) + * + * dB_{0,d}/dc = -e0/2 * T_d(c) + e0 * T_d'(c) for d >= 2 + * + * where T_d(c) is the d >= 2 sum of lcd_binf_b0_zero() without the e0 factor + * and T_d'(c) is its termwise derivative, i.e. term_m scaled by + * -(d - m + 1) / c. Both sums share the one downward recurrence. + * + * @param d index of the base case + * @param e0 exp(-c / 2) + * @param ei0 Ei(-c / 2) + * @param c squared norm of the sample, must be > 0 + * @return d/dc B_{0,d}(0, c) + */ +inline double lcd_binf_b0_zero_deriv(size_t d, double e0, double ei0, + double c) { + assert(c > 0.00); + + if (d == 0) return 0.125 * ei0; + if (d == 1) return -0.25 * e0 / c; + + double term = 0.50 / c; // m = d + double sum = term; // T_d(c) + double sumDeriv = -term / c; // T_d'(c), factor -(d - m + 1) = -1 at m = d + for (size_t m = d; m > 2; --m) { + term *= 2.00 * static_cast(m - 2) / c; // step down to m - 1 + sum += term; + sumDeriv -= static_cast(d - m + 2) * term / c; + } + + return e0 * (sumDeriv - 0.50 * sum); +} + +/** + * @brief B_{k,k}(0, c), the lower-limit attraction term + * + * B_{k,k}(0,c) = 2^-k * sum_{j=0..k} (-1)^j * binom(k,j) * B_{0,j}(0,c) + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param c squared norm of the sample, must be > 0 + * @return B_{k,k}(0, c) + * @note standard-normal target and even N only + */ +inline double lcd_binf_bkk_zero(size_t k, double c) { + assert(k >= 1); + assert(c > 0.00); + + const double e0 = std::exp(-0.50 * c); + const double ei0 = lcd_ei(-0.50 * c); + + double sum = 0.00; + double binomial = 1.00; // binom(k, j) + for (size_t j = 0; j <= k; ++j) { + const double term = lcd_binf_b0_zero(j, e0, ei0, 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 d/dc B_{k,k}(0, c) + * + * dB_{k,k}(0,c)/dc = 2^-k * sum_{j=0..k} (-1)^j * binom(k,j) + * * dB_{0,j}(0,c)/dc + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param c squared norm of the sample, must be > 0 + * @return d/dc B_{k,k}(0, c) + * @note standard-normal target and even N only + */ +inline double lcd_binf_bkk_zero_deriv(size_t k, double c) { + assert(k >= 1); + assert(c > 0.00); + + const double e0 = std::exp(-0.50 * c); + const double ei0 = lcd_ei(-0.50 * c); + + double sum = 0.00; + double binomial = 1.00; // binom(k, j) + for (size_t j = 0; j <= k; ++j) { + const double term = lcd_binf_b0_zero_deriv(j, e0, ei0, 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 S_k(c), the pole-cancelling companion of B_{k,k}(0, c) + * + * S_k(c) = sum_{j=2..k} (-1)^j * binom(k,j) * (j-2)! * 2^(j-3) / c^(j-1) + * + * The sum is empty for k = 1. binom(k,j), (j-2)! * 2^(j-3) and c^(j-1) are + * each carried forward by their own recurrence, as lcd_delta_bkk_zero() does. + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param c squared norm of the sample, must be > 0 + * @return S_k(c), 0 for k = 1 + */ +inline double lcd_binf_s(size_t k, double c) { + assert(k >= 1); + assert(c > 0.00); + + double sum = 0.00; + double binomial = static_cast(k); // binom(k, 1), stepped up below + double coefficient = 0.50; // (j-2)! * 2^(j-3), j = 2 + double cPower = c; // c^(j-1), j = 2 + + for (size_t j = 2; j <= k; ++j) { + binomial *= static_cast(k - j + 1) / static_cast(j); + const double term = binomial * coefficient / cPower; + sum += (j % 2 == 0) ? term : -term; + + coefficient *= 2.00 * static_cast(j - 1); + cPower *= c; + } + + return sum; +} + +/** + * @brief S_k'(c) + * + * S_k'(c) = sum_{j=2..k} (-1)^j * binom(k,j) * (j-2)! * 2^(j-3) + * * (1 - j) / c^j + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param c squared norm of the sample, must be > 0 + * @return S_k'(c), 0 for k = 1 + */ +inline double lcd_binf_s_deriv(size_t k, double c) { + assert(k >= 1); + assert(c > 0.00); + + double sum = 0.00; + double binomial = static_cast(k); // binom(k, 1), stepped up below + double coefficient = 0.50; // (j-2)! * 2^(j-3), j = 2 + double cPower = c; // c^(j-1), j = 2 + + for (size_t j = 2; j <= k; ++j) { + binomial *= static_cast(k - j + 1) / static_cast(j); + const double term = binomial * coefficient * + (1.00 - static_cast(j)) / (cPower * c); + sum += (j % 2 == 0) ? term : -term; + + coefficient *= 2.00 * static_cast(j - 1); + cPower *= c; + } + + return sum; +} + +/** + * @brief phi_k(c), the per-sample attraction potential of the bMax-free + * distance + * + * phi_k(c) = (c + 2k)/8 * ln(c/4) + S_k(c) - 2^k * B_{k,k}(0,c) + * + * The distance uses it as D = K_k - 2 * sum_i w_i * phi_k(c_i) + repulsion. + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param c squared norm of the sample, must be > 0 + * @return phi_k(c) + */ +inline double lcd_binf_phi(size_t k, double c) { + assert(k >= 1); + assert(c > 0.00); + + return 0.125 * (c + 2.00 * static_cast(k)) * std::log(0.25 * c) + + lcd_binf_s(k, c) - + std::ldexp(lcd_binf_bkk_zero(k, c), static_cast(k)); +} + +/** + * @brief phi_k'(c) + * + * phi_k'(c) = ln(c/4)/8 + (c + 2k)/(8c) + S_k'(c) + * - 2^k * dB_{k,k}(0,c)/dc + * + * @param k half the dimension, N = 2k; must be >= 1 + * @param c squared norm of the sample, must be > 0 + * @return phi_k'(c) + */ +inline double lcd_binf_phi_deriv(size_t k, double c) { + assert(k >= 1); + assert(c > 0.00); + + return 0.125 * std::log(0.25 * c) + + (c + 2.00 * static_cast(k)) / (8.00 * c) + + lcd_binf_s_deriv(k, c) - + std::ldexp(lcd_binf_bkk_zero_deriv(k, c), static_cast(k)); +} + +/** + * @brief K_k, the sample-independent constant of the bMax-free distance + * + * K_k = k * sum_{i=2..k} 1/(2i) - 1/2 - k * gamma / 2 + * + * K_k does not depend on the samples, so it shifts the objective without + * moving the optimum. It is carried anyway so the reported value is the true + * distance. + * + * @param k half the dimension, N = 2k; must be >= 1 + * @return K_k + */ +inline double lcd_binf_constant(size_t k) { + assert(k >= 1); + + double harmonic = 0.00; // sum_{i=2..k} 1 / (2i) + for (size_t i = 2; i <= k; ++i) harmonic += 0.50 / static_cast(i); + + return static_cast(k) * harmonic - 0.50 - + 0.50 * static_cast(k) * lcdBinfEulerGamma; +} + +#endif // LCD_EVEN_BINF_CLOSED_FORM_H diff --git a/plots/CMakeLists.txt b/plots/CMakeLists.txt index cecf92e..ed9a3b7 100644 --- a/plots/CMakeLists.txt +++ b/plots/CMakeLists.txt @@ -6,4 +6,8 @@ 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 +target_link_libraries(generate_samples_even PRIVATE approxLCD_static) + +add_executable(generate_samples_even_binf generate_samples_even_binf.cpp) + +target_link_libraries(generate_samples_even_binf PRIVATE approxLCD_static) \ No newline at end of file diff --git a/plots/generate_samples_even_binf.cpp b/plots/generate_samples_even_binf.cpp new file mode 100644 index 0000000..91e5bb1 --- /dev/null +++ b/plots/generate_samples_even_binf.cpp @@ -0,0 +1,73 @@ +#include +#include +#include +#include +#include + +#include "gm_to_dirac_even_binf.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 != 4) { + std::cerr + << "usage: " << argv[0] << " \n" + << " N dimension, must be even and >= 2\n" + << " L number of Dirac components, must be >= 2\n" + << "\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]; + + 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 < 2) { + std::cerr << "error: L must be >= 2; L = 1 has no non-degenerate " + "zero-mean sample set\n"; + return 1; + } + + std::vector x(L * N, 0.00); + GslminimizerResult result; + gm_to_dirac_even_binf approx; + + if (!approx.approximate(L, N, x.data(), nullptr, &result, {})) { + std::cerr << "error: LCD approximation did not converge for N=" << N + << ", L=" << L << "\n"; + return 1; + } + + double distance = 0.00; + approx.modified_van_mises_distance_sq(&distance, L, N, x.data(), nullptr); + + 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 << "distance = " << std::setprecision(17) << distance << '\n'; + std::cerr << "converged in " << result.iterations << " iterations (" + << result.elapsedTimeMicro << " us)\n"; + return 0; +} diff --git a/plots/plot_samples_even_binf.py b/plots/plot_samples_even_binf.py new file mode 100644 index 0000000..7bb95ce --- /dev/null +++ b/plots/plot_samples_even_binf.py @@ -0,0 +1,136 @@ +import argparse +import os +import re +import subprocess +import sys +import tempfile + +import numpy as np + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +GRID_RES = 400 +SEED = "42" + +STANDARD_NORMAL_CASES = [ + {"L": 100}, + {"L": 200}, + {"L": 400}, +] + +DISTANCE_RE = re.compile(r"distance\s*=\s*([-+0-9.eE]+)") + + +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_binf_generator(generator, n_samples, n_dim=2): + """Invoke the bMax-free sampler.""" + 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] + done = subprocess.run( + cmd, env=env, check=True, stderr=subprocess.PIPE, text=True + ) + samples = np.loadtxt(out_csv, delimiter=",") + + match = DISTANCE_RE.search(done.stderr or "") + distance = float(match.group(1)) if match else float("nan") + return samples.reshape(n_samples, n_dim), distance + + +def gaussian_density(axis_limit): + axis = np.linspace(-axis_limit, axis_limit, GRID_RES) + xx, yy = np.meshgrid(axis, axis) + return np.exp(-0.5 * (xx * xx + yy * yy)) / (2.0 * np.pi) + + +def draw_panel(ax, points, title, axis_limit=4.0): + ax.imshow( + gaussian_density(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=10.5) + + +def build_standard_normal_figure(generator, out_path): + fig, axes = plt.subplots( + 1, len(STANDARD_NORMAL_CASES), figsize=(4.2 * len(STANDARD_NORMAL_CASES), 5.1) + ) + fig.patch.set_facecolor("white") + + mean_norms = [] + for ax, case in zip(np.atleast_1d(axes), STANDARD_NORMAL_CASES): + points, distance = run_binf_generator(generator, case["L"]) + covariance = np.cov(points.T, bias=True) + mean_norms.append(np.linalg.norm(points.mean(axis=0))) + draw_panel( + ax, + points, + "L = %d\nΣ = I\nD = %.3e, tr(Cov)/2 = %.4f" + % (case["L"], distance, np.trace(covariance) / 2.0), + ) + + fig.suptitle( + "bMax-free closed-form LCD (N = 2): standard normal with increasing " + "sample count", + fontsize=13, + y=0.985, + ) + fig.tight_layout(rect=(0, 0, 1, 0.90)) + 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_binf"), + help="path to the compiled generate_samples_even_binf binary", + ) + parser.add_argument( + "--out-dir", + default=os.path.join("doxygen", "images"), + help="directory to write the PNG figure into", + ) + args = parser.parse_args() + + generator = resolve_generator(args.generator) + os.makedirs(args.out_dir, exist_ok=True) + + build_standard_normal_figure( + generator, + os.path.join(args.out_dir, "samples_even_binf_standard_normal.png"), + ) + + +if __name__ == "__main__": + main()