Add Von Mises perturbation type for angular parameters

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 &current_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;
 	}

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
 };

source/src/numeric/random/random.cc

@@ -24,6 +24,9 @@
 // C++ 11 Compatibility
 #include <utility/thread/backwards_thread_local.hh>
 
+// Numeric headers
+#include <numeric/constants.hh>
+
 // C++ headers
 #include <iostream>
 #include <cmath>
@@ -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;
 

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