From 8a4d4521758f23c960986a0f59d5925eb1a4292b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 5 Sep 2018 21:09:46 -0400 Subject: made the perturbation method for orthogonal transformations more sensible --- lib/orthogonal.h | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 85d0e11..829df10 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -163,22 +163,36 @@ orthogonal_t generate_rotation_perturbation (gsl_rng *r, const vecto double m_dot_v = 0; for (q_t i = 0; i < q; i++) { - m[0][i] = gsl_ran_ugaussian(r); + m[0][i] = gsl_ran_ugaussian(r); // create a random vector m_dot_v += m[0][i] * v[i]; } double v2 = 0; - double factor = epsilon * gsl_ran_ugaussian(r); for (q_t i = 0; i < q; i++) { - m[0][i] = m[0][i] - m_dot_v * v[i] + factor * v[i]; + m[0][i] = m[0][i] - m_dot_v * v[i]; // find the component orthogonal to v v2 += pow(m[0][i], 2); } double mag_v = sqrt(v2); for (q_t i = 0; i < q; i++) { - m[0][i] /= mag_v; + m[0][i] /= mag_v; // normalize + } + + v2 = 0; + + double factor = epsilon * gsl_ran_ugaussian(r); + + for (q_t i = 0; i < q; i++) { + m[0][i] = m[0][i] + factor * v[i]; // perturb orthogonal vector in original direction + v2 += pow(m[0][i], 2); + } + + mag_v = sqrt(v2); + + for (q_t i = 0; i < q; i++) { + m[0][i] /= mag_v; // normalize } return m; -- cgit v1.2.3-70-g09d2