diff options
| -rw-r--r-- | src/analyze_correlations.cpp | 98 |
1 files changed, 86 insertions, 12 deletions
diff --git a/src/analyze_correlations.cpp b/src/analyze_correlations.cpp index 91d35cb..e108f98 100644 --- a/src/analyze_correlations.cpp +++ b/src/analyze_correlations.cpp @@ -16,6 +16,24 @@ double mean(int N, T *data) { return total / N; } +double squared_mean(int N, double *data) { + double total = 0; + for (int i = 0; i < N; i++) { + total += data[i]; + } + + return total / N; +} + +double central_moment(int N, double *data, double mean, int m) { + double total = 0; + for (int i = 0; i < N; i++) { + total += pow(data[i] - mean, m); + } + + return total / N; +} + void compute_OO(int N, fftw_plan forward_plan, double *forward_data, fftw_plan reverse_plan, double *reverse_data) { fftw_execute(forward_plan); @@ -163,6 +181,10 @@ int main (int argc, char *argv[]) { if (speedy_drop) { drop = N - pow(2, floor(log(N) / log(2))); + } else { + if (N % 2 == 1 && drop % 2 == 0) { + drop++; // make sure M is even + } } if (N <= drop) { @@ -170,15 +192,41 @@ int main (int argc, char *argv[]) { } else { int M = N - drop; + double M_f = (double)M; + + if (length > M) { + length = M; + } + + double *forward_data = (double *)fftw_malloc(M * sizeof(double)); + fftw_plan forward_plan = fftw_plan_r2r_1d(M, forward_data, forward_data, FFTW_R2HC, 0); + double *reverse_data = (double *)fftw_malloc(M * sizeof(double)); + fftw_plan reverse_plan = fftw_plan_r2r_1d(M, reverse_data, reverse_data, FFTW_HC2R, 0); + + uint32_t *data_S = (uint32_t *)malloc(N * sizeof(uint32_t)); fread(data_S, N, sizeof(uint32_t), file_S); - double mean_S = mean(M, data_S + drop); + for (count_t i = 0; i < M; i++) { + forward_data[i] = (double)data_S[drop + i]; + } free(data_S); + double mean_S = mean(M, forward_data); + double squaredMean_S = squared_mean(M, forward_data); + double moment2_S = central_moment(M, forward_data, mean_S, 2); + double moment4_S = central_moment(M, forward_data, mean_S, 4); - double *forward_data = (double *)fftw_malloc(M * sizeof(double)); - fftw_plan forward_plan = fftw_plan_r2r_1d(M, forward_data, forward_data, FFTW_R2HC, 0); + compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); - double M_f = (double)M; + sprintf(filename_S, "wolff_%lu_S_OO.dat", id); + + FILE *file_S = fopen(filename_S, "wb"); + fwrite(&M_f, sizeof(double), 1, file_S); + fwrite(&mean_S, sizeof(double), 1, file_S); + fwrite(&squaredMean_S, sizeof(double), 1, file_S); + fwrite(&moment2_S, sizeof(double), 1, file_S); + fwrite(&moment4_S, sizeof(double), 1, file_S); + fwrite(reverse_data, sizeof(double), length, file_S); + fclose(file_S); uint32_t *data_B = (uint32_t *)malloc((nb - 1) * N * sizeof(uint32_t)); uint32_t *data_M = (uint32_t *)malloc((q - 1) * N * sizeof(uint32_t)); @@ -190,25 +238,23 @@ int main (int argc, char *argv[]) { } double mean_E = mean(M, forward_data); + double squaredMean_E = squared_mean(M, forward_data); + double moment2_E = central_moment(M, forward_data, mean_E, 2); + double moment4_E = central_moment(M, forward_data, mean_E, 4); free(data_B); free(data_M); - double *reverse_data = (double *)fftw_malloc(M * sizeof(double)); - fftw_plan reverse_plan = fftw_plan_r2r_1d(M, reverse_data, reverse_data, FFTW_HC2R, 0); - compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); sprintf(filename_B, "wolff_%lu_E_OO.dat", id); - if (length > M) { - length = M; - } - FILE *file_E = fopen(filename_B, "wb"); fwrite(&M_f, sizeof(double), 1, file_E); fwrite(&mean_E, sizeof(double), 1, file_E); - fwrite(&mean_S, sizeof(double), 1, file_E); + fwrite(&squaredMean_E, sizeof(double), 1, file_E); + fwrite(&moment2_E, sizeof(double), 1, file_E); + fwrite(&moment4_E, sizeof(double), 1, file_E); fwrite(reverse_data, sizeof(double), length, file_E); fclose(file_E); @@ -258,6 +304,10 @@ int main (int argc, char *argv[]) { if (speedy_drop) { drop = N - pow(2, floor(log(N) / log(2))); + } else { + if (N % 2 == 1 && drop % 2 == 0) { + drop++; // make sure M is even + } } if (N <= drop) { @@ -288,6 +338,9 @@ int main (int argc, char *argv[]) { free(data_S); double mean_S = mean(M, forward_data); + double squaredMean_S = squared_mean(M, forward_data); + double moment2_S = central_moment(M, forward_data, mean_S, 2); + double moment4_S = central_moment(M, forward_data, mean_S, 4); compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); @@ -295,6 +348,9 @@ int main (int argc, char *argv[]) { FILE *file_S_new = fopen(filename_S, "wb"); fwrite(&M_f, sizeof(double), 1, file_S_new); fwrite(&mean_S, sizeof(double), 1, file_S_new); + fwrite(&squaredMean_S, sizeof(double), 1, file_S_new); + fwrite(&moment2_S, sizeof(double), 1, file_S_new); + fwrite(&moment4_S, sizeof(double), 1, file_S_new); fwrite(reverse_data, sizeof(double), length, file_S_new); fclose(file_S_new); } @@ -310,6 +366,9 @@ int main (int argc, char *argv[]) { free(data_F); double mean_F = mean(M, forward_data); + double squaredMean_F = squared_mean(M, forward_data); + double moment2_F = central_moment(M, forward_data, mean_F, 2); + double moment4_F = central_moment(M, forward_data, mean_F, 4); compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); @@ -317,6 +376,9 @@ int main (int argc, char *argv[]) { FILE *file_F_new = fopen(filename_F, "wb"); fwrite(&M_f, sizeof(double), 1, file_F_new); fwrite(&mean_F, sizeof(double), 1, file_F_new); + fwrite(&squaredMean_F, sizeof(double), 1, file_F_new); + fwrite(&moment2_F, sizeof(double), 1, file_F_new); + fwrite(&moment4_F, sizeof(double), 1, file_F_new); fwrite(reverse_data, sizeof(double), length, file_F_new); fclose(file_F_new); } @@ -332,6 +394,9 @@ int main (int argc, char *argv[]) { free(data_E); double mean_E = mean(M, forward_data); + double squaredMean_E = squared_mean(M, forward_data); + double moment2_E = central_moment(M, forward_data, mean_E, 2); + double moment4_E = central_moment(M, forward_data, mean_E, 4); compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); @@ -339,6 +404,9 @@ int main (int argc, char *argv[]) { FILE *file_E_new = fopen(filename_E, "wb"); fwrite(&M_f, sizeof(double), 1, file_E_new); fwrite(&mean_E, sizeof(double), 1, file_E_new); + fwrite(&squaredMean_E, sizeof(double), 1, file_E_new); + fwrite(&moment2_E, sizeof(double), 1, file_E_new); + fwrite(&moment4_E, sizeof(double), 1, file_E_new); fwrite(reverse_data, sizeof(double), length, file_E_new); fclose(file_E_new); } @@ -365,6 +433,9 @@ int main (int argc, char *argv[]) { } double mean_M = mean(M, forward_data); + double squaredMean_M = squared_mean(M, forward_data); + double moment2_M = central_moment(M, forward_data, mean_M, 2); + double moment4_M = central_moment(M, forward_data, mean_M, 4); compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); @@ -372,6 +443,9 @@ int main (int argc, char *argv[]) { FILE *file_M_new = fopen(filename_M, "wb"); fwrite(&M_f, sizeof(double), 1, file_M_new); fwrite(&mean_M, sizeof(double), 1, file_M_new); + fwrite(&squaredMean_M, sizeof(double), 1, file_M_new); + fwrite(&moment2_M, sizeof(double), 1, file_M_new); + fwrite(&moment4_M, sizeof(double), 1, file_M_new); fwrite(reverse_data, sizeof(double), length, file_M_new); fclose(file_M_new); } |