From 6ad69f15de3875b6d4f8856885c5064552f64687 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Thu, 30 Apr 2026 15:56:28 +0000 Subject: [PATCH] Add Von Mises perturbation type for angular parameters --- .../parametric/RealValuedParameter.cc | 15 ++++++-- .../parametric/RealValuedParameter.hh | 1 + source/src/numeric/random/random.cc | 37 +++++++++++++++++++ source/src/numeric/random/random.hh | 8 ++++ 4 files changed, 57 insertions(+), 4 deletions(-) diff --git a/source/src/core/conformation/parametric/RealValuedParameter.cc b/source/src/core/conformation/parametric/RealValuedParameter.cc index c12831df931..0aaad7d4e5b 100644 --- a/source/src/core/conformation/parametric/RealValuedParameter.cc +++ b/source/src/core/conformation/parametric/RealValuedParameter.cc @@ -351,8 +351,8 @@ RealValuedParameter::provide_xsd_perturbation_information( ) const { using namespace utility::tag; xsd + XMLSchemaAttribute::attribute_w_default( parameter_name() + "_perturbation" , xsct_real, "Perturbation magnitude for perturbing " + parameter_name() + ".", "0.0" ); - xsd + XMLSchemaAttribute::attribute_w_default( parameter_name() + "_perturbation_type" , xs_string, "Perturbation type for perturbing " + parameter_name() + ". Can be \"gaussian\" or \"uniform\".", "gaussian" ); - //TODO: FIGURE OUT HOW TO ADD RESTRICTION THAT THE PERTURBATION TYPE CAN ONLY BE "gaussian" OR "uniform" + xsd + XMLSchemaAttribute::attribute_w_default( parameter_name() + "_perturbation_type" , xs_string, "Perturbation type for perturbing " + parameter_name() + ". Can be \"gaussian\", \"uniform\", or \"von_mises\". For von_mises, the perturbation magnitude is interpreted as the concentration parameter kappa (higher kappa = tighter distribution; kappa=0 gives uniform on the circle).", "gaussian" ); + //TODO: FIGURE OUT HOW TO ADD RESTRICTION THAT THE PERTURBATION TYPE CAN ONLY BE "gaussian", "uniform", OR "von_mises" } /// @brief Return the XSD information for setting this parameter. @@ -383,7 +383,8 @@ RealValuedParameter::perturbation_type_string_from_enum( ) { static utility::vector1< std::string > const types { //Must match enum in RealValuedParameter.hh. "gaussian", - "uniform" + "uniform", + "von_mises" }; runtime_assert(type_enum > 0 && type_enum < RPT_end_of_list); return types[type_enum]; @@ -406,10 +407,16 @@ RealValuedParameter::perturbation_type_enum_from_string( /// @note Does *not* alter the value of this parameter. core::Real RealValuedParameter::generate_perturbed_value( core::Real const ¤t_value ) const { - runtime_assert_string_msg( perturbation_type_ == RPT_uniform || perturbation_type_ == RPT_gaussian, "Error in RealValuedParameter::generate_perturbed_value(): An unknown perturbation type was set." ); + runtime_assert_string_msg( perturbation_type_ == RPT_uniform || perturbation_type_ == RPT_gaussian || perturbation_type_ == RPT_von_mises, + "Error in RealValuedParameter::generate_perturbed_value(): An unknown perturbation type was set." ); core::Real pert; if ( perturbation_type_ == RPT_uniform ) { pert = (numeric::random::uniform() * 2.0 - 1.0) * perturbation_magnitude_ + current_value; + } else if ( perturbation_type_ == RPT_von_mises ) { + // Von Mises perturbation: perturbation_magnitude_ is the concentration parameter kappa. + // numeric::random::von_mises(kappa) returns a sample centered at 0. + pert = numeric::random::von_mises( perturbation_magnitude_ ) + current_value; + pert = numeric::principal_angle_radians( pert ); } else { // RPT_gaussian case pert = numeric::random::gaussian() * perturbation_magnitude_ + current_value; } diff --git a/source/src/core/conformation/parametric/RealValuedParameter.hh b/source/src/core/conformation/parametric/RealValuedParameter.hh index 7dfd10fd18b..5586180b79c 100644 --- a/source/src/core/conformation/parametric/RealValuedParameter.hh +++ b/source/src/core/conformation/parametric/RealValuedParameter.hh @@ -49,6 +49,7 @@ namespace parametric { enum RealPerturbationType { RPT_gaussian=1, //Keep first RPT_uniform, + RPT_von_mises, RPT_unknown_type, //Keep second-to-last RPT_end_of_list=RPT_unknown_type //Keep last }; diff --git a/source/src/numeric/random/random.cc b/source/src/numeric/random/random.cc index 2353453b546..e2f0456e349 100644 --- a/source/src/numeric/random/random.cc +++ b/source/src/numeric/random/random.cc @@ -24,6 +24,9 @@ // C++ 11 Compatibility #include +// Numeric headers +#include + // C++ headers #include #include @@ -52,6 +55,40 @@ double uniform() { return rg().uniform(); } double gaussian() { return rg().gaussian(); } int random_range(int low, int high) { return rg().random_range(low, high); } +/// @brief Generate a random deviate from the Von Mises distribution centered at 0 +/// with concentration parameter kappa, using the Best-Fisher (1979) algorithm. +/// @details When kappa is 0, returns a uniform random angle in (-pi, pi]. +/// For large kappa, the distribution approximates a Gaussian with sigma = 1/sqrt(kappa). +/// @author Andy Watkins +double von_mises( double kappa ) { + double const pi = numeric::constants::d::pi; + + // Special case: kappa == 0 gives uniform distribution on the circle + if ( kappa < 1e-12 ) { + return ( uniform() * 2.0 - 1.0 ) * pi; + } + + // Best-Fisher algorithm (Best & Fisher, 1979, Applied Statistics 28:152-157) + double const tau = 1.0 + std::sqrt( 1.0 + 4.0 * kappa * kappa ); + double const rho = ( tau - std::sqrt( 2.0 * tau ) ) / ( 2.0 * kappa ); + double const r = ( 1.0 + rho * rho ) / ( 2.0 * rho ); + + while ( true ) { + double const u1 = uniform(); + double const z = std::cos( pi * u1 ); + double const f = ( 1.0 + r * z ) / ( r + z ); + double const c = kappa * ( r - f ); + + double const u2 = uniform(); + + if ( c * ( 2.0 - c ) > u2 || std::log( c / u2 ) + 1.0 >= c ) { + double const u3 = uniform(); + double const theta = ( u3 > 0.5 ) ? std::acos( f ) : -std::acos( f ); + return theta; + } + } +} + using namespace std; diff --git a/source/src/numeric/random/random.hh b/source/src/numeric/random/random.hh index 0c8709597d5..8a97411aa45 100644 --- a/source/src/numeric/random/random.hh +++ b/source/src/numeric/random/random.hh @@ -68,6 +68,14 @@ double gaussian(); /// and high. Threadsafe since each thread uses its own random generator. int random_range(int low, int high); +/// @brief Generate a random number drawn from the Von Mises distribution +/// centered at 0 with concentration parameter kappa. The Von Mises distribution +/// is the circular analog of the Gaussian distribution; for large kappa it +/// approximates a Gaussian with sigma = 1/sqrt(kappa). When kappa is 0, this +/// returns a uniform random angle in (-pi, pi]. Uses the Best-Fisher (1979) +/// algorithm. Threadsafe since each thread uses its own random generator. +double von_mises( double kappa ); + /// @brief Random number generator system class RandomGenerator : public utility::VirtualBase