diff options
| -rw-r--r-- | src/analyze_correlations.cpp | 95 |
1 files changed, 71 insertions, 24 deletions
diff --git a/src/analyze_correlations.cpp b/src/analyze_correlations.cpp index e5cd200..2170a70 100644 --- a/src/analyze_correlations.cpp +++ b/src/analyze_correlations.cpp @@ -6,50 +6,51 @@ #include <getopt.h> #include <fftw3.h> -double *compute_OO(count_t length, double *tmp_data, count_t N) { +double *compute_OO(count_t length, double *data, count_t N) { double *OO = (double *)calloc(2 + length, sizeof(double)); if (OO == NULL) { - printf("calloc failure\n"); - exit(EXIT_FAILURE); + return NULL; } -// int better_N = pow(2, floor(log(N) / log(2))); - int better_N = N; - double *data = tmp_data + (N - better_N); - - OO[0] = (double)better_N; + OO[0] = (double)N; - for (count_t i = 0; i < better_N; i++) { + for (count_t i = 0; i < N; i++) { OO[1] += data[i]; } - OO[1] /= better_N; + OO[1] /= N; - fftw_plan plan = fftw_plan_r2r_1d(better_N, data, data, FFTW_R2HC, FFTW_ESTIMATE); + fftw_plan plan = fftw_plan_r2r_1d(N, data, data, FFTW_R2HC, FFTW_ESTIMATE); fftw_execute(plan); - double *tmp = (double *)malloc(better_N * sizeof(double)); + double *tmp = (double *)malloc(N * sizeof(double)); + + if (tmp == NULL) { + free(OO); + fftw_destroy_plan(plan); + return NULL; + } tmp[0] = data[0] * data[0]; - tmp[better_N / 2] = data[better_N/2] * data[better_N/2]; + tmp[N / 2] = data[N/2] * data[N/2]; - for (count_t i = 1; i < better_N / 2 - 1; i++) { - tmp[i] = data[i] * data[i] + data[better_N - i] * data[better_N - i]; - tmp[better_N - i] = 0; + for (count_t i = 1; i < N / 2 - 1; i++) { + tmp[i] = data[i] * data[i] + data[N - i] * data[N - i]; + tmp[N - i] = 0; } fftw_destroy_plan(plan); - fftw_plan plan2 = fftw_plan_r2r_1d(better_N, tmp, tmp, FFTW_HC2R, FFTW_ESTIMATE); - fftw_execute(plan2); + plan = fftw_plan_r2r_1d(N, tmp, tmp, FFTW_HC2R, FFTW_ESTIMATE); + fftw_execute(plan); for (count_t j = 0; j < length; j++) { - OO[2 + j] = tmp[j] / pow(better_N, 2); + OO[2 + j] = tmp[j] / pow(N, 2); } free(tmp); - fftw_destroy_plan(plan2); + fftw_destroy_plan(plan); return OO; } @@ -79,10 +80,12 @@ double finite_energy(q_t nb, double *J, q_t q, double *H, v_t nv, v_t ne, uint32 int main (int argc, char *argv[]) { count_t drop = (count_t)1e4; count_t length = (count_t)1e4; + bool speedy_drop = false; + bool from_stdin = false; int opt; - while ((opt = getopt(argc, argv, "d:l:")) != -1) { + while ((opt = getopt(argc, argv, "d:l:sp")) != -1) { switch (opt) { case 'd': drop = (count_t)atof(optarg); @@ -90,19 +93,37 @@ int main (int argc, char *argv[]) { case 'l': length = (count_t)atof(optarg); break; + case 's': + speedy_drop = true; + break; + case 'p': + from_stdin = true; + break; default: exit(EXIT_FAILURE); } } + FILE *metadata; - FILE *metadata = fopen("wolff_metadata.txt", "r"); + if (from_stdin) { + metadata = stdin; + } else { + metadata = fopen("wolff_metadata.txt", "r"); + } if (metadata == NULL) { + printf("Metadata file not found. Make sure you are in the correct directory!\n"); exit(EXIT_FAILURE); } unsigned long id; char *model = (char *)malloc(32 * sizeof(char)); + + if (model == NULL) { + printf("Malloc failed.\n"); + exit(EXIT_FAILURE); + } + q_t q; D_t D; L_t L; @@ -110,13 +131,20 @@ int main (int argc, char *argv[]) { while (EOF != fscanf(metadata, "<| \"ID\" -> %lu, \"MODEL\" -> \"%[^\"]\", \"q\" -> %" SCNq ", \"D\" -> %" SCND ", \"L\" -> %" SCNL ", \"NV\" -> %" SCNv ", \"NE\" -> %" SCNv ", ", &id, model, &q, &D, &L, &nv, &ne)) { - if (0 == strcmp(model, "ISING") || 0 == strcmp(model, "POTTS") || 0 == strcmp(model, "CLOCK") || 0 == strcmp(model, "DGM")) { + bool is_finite = 0 == strcmp(model, "ISING") || 0 == strcmp(model, "POTTS") || 0 == strcmp(model, "CLOCK") || 0 == strcmp(model, "DGM"); + + if (is_finite) { q_t nb; double T; fscanf(metadata, "\"NB\" -> %" SCNq ", \"T\" -> %lf, \"J\" -> {", &nb, &T); double *J = (double *)malloc(nb * sizeof(double)); double *H = (double *)malloc(q * sizeof(double)); + if (J == NULL || H == NULL) { + printf("%lu: Malloc failed.\n", id); + break; + } + for (q_t i = 0; i < nb - 1; i++) { fscanf(metadata, "%lf, ", &(J[i])); } @@ -130,6 +158,11 @@ int main (int argc, char *argv[]) { char *filename_B = (char *)malloc(128 * sizeof(char)); char *filename_S = (char *)malloc(128 * sizeof(char)); + if (filename_M == NULL || filename_B == NULL || filename_S == NULL) { + printf("%lu: Malloc failed.\n", id); + break; + } + sprintf(filename_M, "wolff_%lu_M.dat", id); sprintf(filename_B, "wolff_%lu_B.dat", id); sprintf(filename_S, "wolff_%lu_S.dat", id); @@ -138,12 +171,22 @@ int main (int argc, char *argv[]) { FILE *file_B = fopen(filename_B, "rb"); FILE *file_S = fopen(filename_S, "rb"); + if (file_M == NULL || file_B == NULL || file_S == NULL) { + printf("%lu: Opening data file failed.\n", id); + break; + } + fseek(file_S, 0, SEEK_END); unsigned long N = ftell(file_S) / sizeof(uint32_t); fseek(file_S, 0, SEEK_SET); + if (speedy_drop) { + drop = N - pow(2, floor(log(N) / log(2))); + } + if (N < drop) { - printf("Number of steps %lu is less than %" PRIcount ", nothing done.\n", N, drop); + printf("%lu: Number of steps %lu is less than %" PRIcount ", nothing done.\n", id, N, drop); + break; } else { double *data_S = (double *)malloc(N * sizeof(double)); @@ -232,6 +275,10 @@ int main (int argc, char *argv[]) { unsigned long N = ftell(file_S) / sizeof(uint32_t); fseek(file_S, 0, SEEK_SET); + if (speedy_drop) { + drop = N - pow(2, floor(log(N) / log(2))); + } + if (N < drop) { printf("Number of steps %lu is less than %" PRIcount ", nothing done.\n", N, drop); } else { |