diff options
| -rw-r--r-- | lib/cluster.c | 91 | ||||
| -rw-r--r-- | lib/cluster_finite.c | 10 | ||||
| -rw-r--r-- | lib/cluster_finite.h | 7 | ||||
| -rw-r--r-- | lib/dihedral.c | 23 | ||||
| -rw-r--r-- | lib/dihedral.h | 4 | ||||
| -rw-r--r-- | lib/initial_finite.c | 145 | ||||
| -rw-r--r-- | lib/symmetric.c | 43 | ||||
| -rw-r--r-- | lib/symmetric.h | 2 | ||||
| -rw-r--r-- | src/wolff_finite.c | 12 |
9 files changed, 207 insertions, 130 deletions
diff --git a/lib/cluster.c b/lib/cluster.c index 7274eb9..96225a2 100644 --- a/lib/cluster.c +++ b/lib/cluster.c @@ -1,97 +1,6 @@ #include "cluster.h" -v_t flip_cluster(ising_state_t *s, v_t v0, q_t rot, gsl_rng *r) { - v_t nv = 0; - - ll_t *stack = NULL; // create a new stack - stack_push(&stack, v0); // push the initial vertex to the stack - - bool *marks = (bool *)calloc(s->g->nv, sizeof(bool)); - - while (stack != NULL) { - v_t v = stack_pop(&stack); - - if (!marks[v]) { - q_t s_old, s_new; - dihedral_t *R_new; - bool external_flipped; - - marks[v] = true; - - if (v == s->g->nv - 1) { - R_new = dihedral_compose(s->q, rot, s->R); - external_flipped = true; - } else { - s_old = s->spins[v]; - s_new = dihedral_act(s->q, rot, s_old); - external_flipped = false; - } - - v_t nn = s->g->v_i[v + 1] - s->g->v_i[v]; - - for (v_t i = 0; i < nn; i++) { - q_t sn; - double prob; - bool external_neighbor = false; - - v_t vn = s->g->v_adj[s->g->v_i[v] + i]; - - if (vn == s->g->nv - 1) { - external_neighbor = true; - } else { - sn = s->spins[vn]; - } - - if (external_flipped || external_neighbor) { - q_t rot_s_old, rot_s_new; - - if (external_neighbor) { - rot_s_old = dihedral_inverse_act(s->q, s->R, s_old); - rot_s_new = dihedral_inverse_act(s->q, s->R, s_new); - } else { - rot_s_old = dihedral_inverse_act(s->q, s->R, sn); - rot_s_new = dihedral_inverse_act(s->q, R_new, sn); - } - - prob = s->H_probs[rot_s_new * s->q + rot_s_old]; - - s->M[rot_s_old]--; - s->M[rot_s_new]++; - - s->E += - s->H[rot_s_new] + s->H[rot_s_old]; - } else { - q_t diff_old = (s_old + s->q - sn) % s->q; - q_t diff_new = (s_new + s->q - sn) % s->q; - - prob = s->J_probs[diff_new * s->q + diff_old]; - - s->E += - s->J[diff_new] + s->J[diff_old]; - } - - if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... - stack_push(&stack, vn); // push the neighboring vertex to the stack - } - } - - if (external_flipped) { - free(s->R); - s->R = R_new; - } else { - s->spins[v] = s_new; - } - - if (v != s->g->nv - 1) { // count the number of non-external sites that flip - nv++; - } - } - } - - free(marks); - - return nv; -} - v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t rot, gsl_rng *r) { v_t nv = 0; diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c index 71396e0..9392cf8 100644 --- a/lib/cluster_finite.c +++ b/lib/cluster_finite.c @@ -1,8 +1,8 @@ #include "cluster_finite.h" -v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) { - q_t *rot = s->transformations + s->q * rot_ind; +v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) { + q_t *rot = s->transformations + s->q * s->involutions[rot_ind]; q_t *R_inv = symmetric_invert(s->q, s->R); v_t nv = 0; @@ -63,10 +63,10 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) { s->M[rot_s_new]++; } else { - q_t diff_old = (s_old + s->q - sn) % s->q; - q_t diff_new = (s_new + s->q - sn) % s->q; + q_t diff_old = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_old]]; + q_t diff_new = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_new]]; - prob = s->J_probs[diff_new * s->q + diff_old]; + prob = s->J_probs[diff_new * s->n_bond_types + diff_old]; s->B[diff_old]--; s->B[diff_new]++; diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h index 701c197..b2d764e 100644 --- a/lib/cluster_finite.h +++ b/lib/cluster_finite.h @@ -31,6 +31,11 @@ typedef struct { q_t q; R_t n_transformations; q_t *transformations; + R_t n_involutions; + R_t *involutions; + R_t *transform_site_to_zero; + q_t n_bond_types; + q_t *bond_with_zero_type; double T; double *J; double *H; @@ -42,5 +47,5 @@ typedef struct { v_t *M; } state_finite_t; -v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot, gsl_rng *r); +v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot, gsl_rng *r); diff --git a/lib/dihedral.c b/lib/dihedral.c index ac74a23..8158b43 100644 --- a/lib/dihedral.c +++ b/lib/dihedral.c @@ -11,10 +11,14 @@ dihedral_t *dihedral_compose(q_t q, q_t g1i, const dihedral_t *g2) { return g3; } -q_t dihedral_act(q_t q, q_t gi, q_t s) { +q_t dihedral_act(q_t q, q_t gi, bool r, q_t s) { // we only need to consider the action of reflections - return (gi + q - s) % q; + if (r) { + return (gi + q - s) % q; + } else { + return (gi + s) % q; + } } q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) { @@ -26,15 +30,26 @@ q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) { } q_t *dihedral_gen_transformations(q_t q) { - q_t *transformations = (q_t *)malloc(q * q * sizeof(q_t)); + q_t *transformations = (q_t *)malloc(2 * q * q * sizeof(q_t)); for (q_t i = 0; i < q; i++) { for (q_t j = 0; j < q; j++) { - transformations[q * i + j] = dihedral_act(q, i, j); + transformations[q * i + j] = dihedral_act(q, i, false, j); + transformations[q * q + q * i + j] = dihedral_act(q, i, true, j); } } return transformations; } +R_t *dihedral_gen_involutions(q_t q) { + R_t *transformations = (R_t *)malloc(q * sizeof(R_t)); + + for (q_t i = 0; i < q; i++) { + transformations[i] = q + i; + } + + return transformations; +} + diff --git a/lib/dihedral.h b/lib/dihedral.h index e5e4cbd..c95b23a 100644 --- a/lib/dihedral.h +++ b/lib/dihedral.h @@ -11,9 +11,11 @@ typedef struct { dihedral_t *dihedral_compose(q_t q, q_t gti, const dihedral_t *g2); -q_t dihedral_act(q_t q, q_t gi, q_t s); +q_t dihedral_act(q_t q, q_t gi, bool r, q_t s); q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s); q_t *dihedral_gen_transformations(q_t q); +R_t *dihedral_gen_involutions(q_t q); +R_t factorial(q_t); diff --git a/lib/initial_finite.c b/lib/initial_finite.c index fb120f0..6ea76ef 100644 --- a/lib/initial_finite.c +++ b/lib/initial_finite.c @@ -23,6 +23,40 @@ q_t *initialize_R(q_t q) { return R; } +R_t *transformation_bringing_to_zero(q_t q, R_t n_transformations, q_t *transformations) { + R_t *destination = (R_t *)malloc(q * sizeof(R_t)); + + for (q_t i = 0; i < q; i++) { + for (R_t j = 0; j < n_transformations; j++) { + if (transformations[q * j + i] == 0) { + destination[i] = j; + } + } + } + + return destination; +} + +R_t find_involutions(R_t *destination, q_t q, R_t n_transformations, q_t *transformations) { + R_t n_involutions = 0; + + for (R_t i = 1; i < n_transformations; i++) { + bool is_involution = true; + for (q_t j = 0; j < q; j++) { + if (j != transformations[q * i + transformations[q * i + j]]) { + is_involution = false; + break; + } + } + if (is_involution) { + destination[n_involutions] = i; + n_involutions++; + } + } + + return n_involutions; +} + state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) { state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t)); @@ -38,11 +72,26 @@ state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) } s->q = 2; - s->n_transformations = 1; - s->transformations = (q_t *)malloc(2 * sizeof(q_t)); - s->transformations[0] = 1; - s->transformations[1] = 0; + s->n_transformations = 2; + s->transformations = (q_t *)malloc(2 * 2 * sizeof(q_t)); + s->transformations[0] = 0; + s->transformations[1] = 1; + s->transformations[2] = 1; + s->transformations[3] = 0; + + s->n_involutions = 1; + s->involutions = (R_t *)malloc(1 * sizeof(R_t)); + s->involutions[0] = 1; + + s->transform_site_to_zero = (R_t *)malloc(2 * sizeof(R_t)); + s->transform_site_to_zero[0] = 0; + s->transform_site_to_zero[1] = 1; + + s->n_bond_types = 2; + s->bond_with_zero_type = (q_t *)malloc(2 * sizeof(q_t)); + s->bond_with_zero_type[0] = 0; + s->bond_with_zero_type[1] = 1; s->T = T; s->J = (double *)malloc(2 * sizeof(double)); @@ -81,19 +130,34 @@ state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, doub } s->q = q; - s->n_transformations = q; - s->transformations = dihedral_gen_transformations(q); + s->n_transformations = factorial(q); + s->transformations = symmetric_gen_transformations(q); + s->involutions = (R_t *)malloc(s->n_transformations * sizeof(R_t)); + s->n_involutions = find_involutions(s->involutions, q, s->n_transformations, s->transformations); + + s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations); + + s->n_bond_types = 2; + + s->bond_with_zero_type = (q_t *)malloc(q * sizeof(q_t)); + + s->bond_with_zero_type[0] = 0; + + for (q_t i = 1; i < q; i++) { + s->bond_with_zero_type[i] = 1; + } s->T = T; - s->J = (double *)calloc(q, sizeof(double)); + s->J = (double *)calloc(2, sizeof(double)); s->J[0] = 1.0; + s->J[1] = 0.0; s->H = (double *)malloc(q * sizeof(double)); for (q_t i = 0; i < q; i++) { s->H[i] = H[i]; } - s->J_probs = Jprobs_from_J(q, T, s->J); + s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J); s->H_probs = Jprobs_from_J(q, T, s->H); s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); @@ -101,7 +165,7 @@ state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, doub s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? - s->B = (v_t *)calloc(q, sizeof(v_t)); + s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); s->B[0] = s->ne; // everyone starts in state 0, remember? return s; @@ -122,13 +186,30 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub } s->q = q; - s->n_transformations = q; + + s->n_transformations = 2 * q; s->transformations = dihedral_gen_transformations(q); + s->n_involutions = q; + s->involutions = dihedral_gen_involutions(q); + + s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations); + s->bond_with_zero_type = malloc(q * sizeof(q_t)); + + s->n_bond_types = q / 2 + 1; + + for (q_t i = 0; i < q / 2 + 1; i++) { + s->bond_with_zero_type[i] = i; + } + + for (q_t i = 1; i < (q + 1) / 2; i++) { + s->bond_with_zero_type[q - i] = i; + } + s->T = T; - s->J = (double *)malloc(q * sizeof(double)); + s->J = (double *)malloc(s->n_bond_types * sizeof(double)); - for (q_t i = 0; i < q; i++) { + for (q_t i = 0; i < s->n_bond_types; i++) { s->J[i] = cos(2 * M_PI * i / ((double)q)); } @@ -138,7 +219,7 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub s->H[i] = H[i]; } - s->J_probs = Jprobs_from_J(q, T, s->J); + s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J); s->H_probs = Jprobs_from_J(q, T, s->H); s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); @@ -146,7 +227,7 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? - s->B = (v_t *)calloc(q, sizeof(v_t)); + s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); s->B[0] = s->ne; // everyone starts in state 0, remember? return s; @@ -168,17 +249,30 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double } s->q = q; - s->n_transformations = q; + + s->n_transformations = 2 * q; s->transformations = dihedral_gen_transformations(q); + s->n_involutions = q; + s->involutions = dihedral_gen_involutions(q); - s->T = T; - s->J = (double *)malloc(q * sizeof(double)); + s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations); + s->bond_with_zero_type = malloc(q * sizeof(q_t)); + + s->n_bond_types = q / 2 + 1; for (q_t i = 0; i < q / 2 + 1; i++) { - s->J[i] = -pow(i, 2); + s->bond_with_zero_type[i] = i; } - for (q_t i = 1; i < (q + 1) / 2; i++) { - s->J[q - i] = -pow(i, 2); + + for (q_t i = 1; i < (q + 1) / 2; i++) { + s->bond_with_zero_type[(int)q - (int)i] = i; + } + + s->T = T; + s->J = (double *)malloc(s->n_bond_types * sizeof(double)); + + for (q_t i = 0; i < s->n_bond_types; i++) { + s->J[i] = -pow(i, 2); } s->H = (double *)malloc(q * sizeof(double)); @@ -186,7 +280,7 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double s->H[i] = H[i]; } - s->J_probs = Jprobs_from_J(q, T, s->J); + s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J); s->H_probs = Jprobs_from_J(q, T, s->H); s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); @@ -194,6 +288,8 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? + s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); + s->B[0] = s->nv; // everyone starts in state 0, remember? return s; } @@ -201,8 +297,10 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double double state_finite_energy(state_finite_t *s) { double E = 0; - for (q_t i = 0; i < s->q; i++) { + for (q_t i = 0; i < s->n_bond_types; i++) { E += s->J[i] * s->B[i]; + } + for (q_t i = 0; i < s->q; i++) { E += s->H[i] * s->M[i]; } @@ -220,6 +318,9 @@ void state_finite_free(state_finite_t *s) { free(s->M); free(s->B); free(s->transformations); + free(s->involutions); + free(s->transform_site_to_zero); + free(s->bond_with_zero_type); free(s); } diff --git a/lib/symmetric.c b/lib/symmetric.c index 729b38c..4487538 100644 --- a/lib/symmetric.c +++ b/lib/symmetric.c @@ -25,3 +25,46 @@ q_t *symmetric_invert(q_t q, const q_t *g) { return g_inv; } +void swap(q_t *q1, q_t *q2) { + q_t temp = *q1; + *q1 = *q2; + *q2 = temp; +} + +R_t factorial(q_t q) { + if (q == 0) { + return 1; + } else { + return q * factorial(q - 1); + } +} + +void permute(q_t *a, q_t l, q_t r, R_t pos, q_t *transformations) { + if (l == r - 1) { + for (q_t i = 0; i < r; i++) { + transformations[r * pos + i] = a[i]; + } + } else { + for (q_t i = l; i < r; i++) { + swap((a+l), (a+i)); + permute(a, l+1, r, pos + (i - l) * factorial(r - l - 1), transformations); + swap((a+l), (a+i)); + } + } +} + +q_t *symmetric_gen_transformations(q_t q) { + q_t *transformations = (q_t *)malloc(q * factorial(q) * sizeof(q_t)); + q_t *tmp = (q_t *)malloc(q * sizeof(q_t)); + + for (q_t i = 0; i < q; i++) { + tmp[i] = i; + } + + permute(tmp, 0, q, 0, transformations); + + free(tmp); + + return transformations; +} + diff --git a/lib/symmetric.h b/lib/symmetric.h index 6e00f52..c71521d 100644 --- a/lib/symmetric.h +++ b/lib/symmetric.h @@ -11,3 +11,5 @@ q_t symmetric_act(const q_t *g, q_t s); q_t *symmetric_invert(q_t q, const q_t *g); +q_t *symmetric_gen_transformations(q_t q); + diff --git a/src/wolff_finite.c b/src/wolff_finite.c index e41c326..4bf96b9 100644 --- a/src/wolff_finite.c +++ b/src/wolff_finite.c @@ -96,7 +96,7 @@ int main(int argc, char *argv[]) { char *filename_S = (char *)malloc(256 * sizeof(char)); sprintf(filename_M, "wolff_%s_%" PRIq "_%" PRID "_%" PRIL "_%.8f", finite_model_t_strings[model], q, D, L, T); - for (q_t i = 0; i < q; i++) { + for (q_t i = 0; i < s->q; i++) { sprintf(filename_M + strlen(filename_M), "_%.8f", s->H[i]); } @@ -126,8 +126,8 @@ int main(int argc, char *argv[]) { bool changed = false; while (!changed) { - step = gsl_rng_uniform_int(r, s->n_transformations); - if (symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) { + step = gsl_rng_uniform_int(r, s->n_involutions); + if (symmetric_act(s->transformations + s->q * s->involutions[step], s->spins[v0]) != s->spins[v0]) { changed = true; } } @@ -138,9 +138,9 @@ int main(int argc, char *argv[]) { // as a fast time many machines will actually have it be 64 bits. we cast // it down here to halve space. - for (q_t i = 0; i < q - 1; i++) { - fwrite(&(s->M[i]), sizeof(uint32_t), 1, outfile_M); // if we know the occupation of the first q - 1 states, we know the occupation of the last - fwrite(&(s->B[i]), sizeof(uint32_t), 1, outfile_B); // if we know the occupation of the first q - 1 states, we know the occupation of the last + for (q_t i = 0; i < s->q - 1; i++) { // if we know the occupation of the first q - 1 states, we know the occupation of the last + fwrite(&(s->M[i]), sizeof(uint32_t), 1, outfile_M); + fwrite(&(s->B[i]), sizeof(uint32_t), 1, outfile_B); } fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); |