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-54-g00ecf From 10510b7e2a09c102c307165e0cd716e17b17ef84 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 5 Sep 2018 22:45:39 -0400 Subject: fixed a bug in perturbative generation and set up epsilon to be automatically set by T and H --- lib/orthogonal.h | 2 +- src/wolff_On.cpp | 16 ++++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 829df10..4ae3ca6 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -182,7 +182,7 @@ orthogonal_t generate_rotation_perturbation (gsl_rng *r, const vecto v2 = 0; - double factor = epsilon * gsl_ran_ugaussian(r); + double factor = gsl_ran_gaussian(r, epsilon); for (q_t i = 0; i < q; i++) { m[0][i] = m[0][i] + factor * v[i]; // perturb orthogonal vector in original direction diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 6cd8777..736eb1b 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -142,6 +142,22 @@ int main(int argc, char *argv[]) { std::function gen_R; if (use_pert) { + double Hish; + if (modulated_field) { + Hish = fabs(H_vec[0]); + } else { + double H2 = 0; + for (q_t i = 0; i < N_COMP; i++) { + H2 += pow(H_vec[i], 2); + } + Hish = sqrt(H2); + } + + if (Hish > 1.0) { + epsilon = sqrt(T / Hish); + } else { + epsilon = sqrt(T); + } gen_R = std::bind(generate_rotation_perturbation , std::placeholders::_1, std::placeholders::_2, epsilon, order); pert_type = "PERTURB"; } else { -- cgit v1.2.3-54-g00ecf From 5d074a4f70dcfc099ddbb5483079b19dca18d103 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 5 Sep 2018 22:48:15 -0400 Subject: right now I need to differentiate between the perturbation types --- src/wolff_On.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 736eb1b..8a7b1a6 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -159,7 +159,7 @@ int main(int argc, char *argv[]) { epsilon = sqrt(T); } gen_R = std::bind(generate_rotation_perturbation , std::placeholders::_1, std::placeholders::_2, epsilon, order); - pert_type = "PERTURB"; + pert_type = "PERTURB2"; } else { gen_R = generate_rotation_uniform ; pert_type = "UNIFORM"; -- cgit v1.2.3-54-g00ecf From baca97666aebc1926861338af312c0560bf923f7 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 6 Sep 2018 23:36:19 -0400 Subject: tweaked the perturbation yet again --- src/wolff_On.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 8a7b1a6..fada9bf 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -153,13 +153,10 @@ int main(int argc, char *argv[]) { Hish = sqrt(H2); } - if (Hish > 1.0) { - epsilon = sqrt(T / Hish); - } else { - epsilon = sqrt(T); - } + epsilon = sqrt(N_COMP * T / (D + Hish / 2)); + gen_R = std::bind(generate_rotation_perturbation , std::placeholders::_1, std::placeholders::_2, epsilon, order); - pert_type = "PERTURB2"; + pert_type = "PERTURB3"; } else { gen_R = generate_rotation_uniform ; pert_type = "UNIFORM"; -- cgit v1.2.3-54-g00ecf From 42efc30ddcb2d0b8471120d56d2c0c64d6866aa9 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 6 Sep 2018 23:41:22 -0400 Subject: small changed to perturb again --- src/wolff_On.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index fada9bf..0de9d42 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -153,10 +153,10 @@ int main(int argc, char *argv[]) { Hish = sqrt(H2); } - epsilon = sqrt(N_COMP * T / (D + Hish / 2)); + epsilon = sqrt((N_COMP - 1) * T / (D + Hish / 2)); gen_R = std::bind(generate_rotation_perturbation , std::placeholders::_1, std::placeholders::_2, epsilon, order); - pert_type = "PERTURB3"; + pert_type = "PERTURB4"; } else { gen_R = generate_rotation_uniform ; pert_type = "UNIFORM"; -- cgit v1.2.3-54-g00ecf From 749d387b48456b78e448ab33516c01de20a2624b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 8 Sep 2018 19:10:08 -0400 Subject: small perturb change --- lib/orthogonal.h | 2 +- src/wolff_On.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 4ae3ca6..e1bf33a 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -185,7 +185,7 @@ orthogonal_t generate_rotation_perturbation (gsl_rng *r, const vecto double factor = gsl_ran_gaussian(r, epsilon); for (q_t i = 0; i < q; i++) { - m[0][i] = m[0][i] + factor * v[i]; // perturb orthogonal vector in original direction + m[0][i] += factor * v[i]; // perturb orthogonal vector in original direction v2 += pow(m[0][i], 2); } diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 0de9d42..2356875 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -153,10 +153,10 @@ int main(int argc, char *argv[]) { Hish = sqrt(H2); } - epsilon = sqrt((N_COMP - 1) * T / (D + Hish / 2)); + epsilon = sqrt((N_COMP - 1) * T / (D + Hish / 2)) / 2; gen_R = std::bind(generate_rotation_perturbation , std::placeholders::_1, std::placeholders::_2, epsilon, order); - pert_type = "PERTURB4"; + pert_type = "PERTURB5"; } else { gen_R = generate_rotation_uniform ; pert_type = "UNIFORM"; -- cgit v1.2.3-54-g00ecf