Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Add and make use of a tiny statistics library for easier handling of …

…the measurements.
  • Loading branch information...
commit 18f6f3b074fb5e2ddd47bfb4657ed30a2338d8c5 1 parent 631aefa
@MrMage MrMage authored
View
6 qed.xcodeproj/project.pbxproj
@@ -7,6 +7,7 @@
objects = {
/* Begin PBXBuildFile section */
+ 6561B0F414867BCF0096B9CB /* statistics.c in Sources */ = {isa = PBXBuildFile; fileRef = 6561B0F314867BCF0096B9CB /* statistics.c */; };
65AC1ABA14698B2B00336FFC /* leapfrog2.c in Sources */ = {isa = PBXBuildFile; fileRef = 65AC1AB814698B2B00336FFC /* leapfrog2.c */; };
65AC1ABE14698B4400336FFC /* rec_lf_integrator.c in Sources */ = {isa = PBXBuildFile; fileRef = 65AC1ABC14698B4400336FFC /* rec_lf_integrator.c */; };
65E7DFF014600AD300095D43 /* gauss.c in Sources */ = {isa = PBXBuildFile; fileRef = 65E7DFD614600AD300095D43 /* gauss.c */; };
@@ -49,6 +50,8 @@
/* End PBXCopyFilesBuildPhase section */
/* Begin PBXFileReference section */
+ 6561B0F314867BCF0096B9CB /* statistics.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = statistics.c; sourceTree = "<group>"; };
+ 6561B0F614867BD70096B9CB /* statistics.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = statistics.h; sourceTree = "<group>"; };
65AC1AB814698B2B00336FFC /* leapfrog2.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = leapfrog2.c; sourceTree = "<group>"; };
65AC1AB914698B2B00336FFC /* leapfrog2.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = leapfrog2.h; sourceTree = "<group>"; };
65AC1ABC14698B4400336FFC /* rec_lf_integrator.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = rec_lf_integrator.c; sourceTree = "<group>"; };
@@ -160,6 +163,8 @@
65E7DFEC14600AD300095D43 /* linalg.h */,
65E7DFED14600AD300095D43 /* measurements.c */,
65E7DFEE14600AD300095D43 /* measurements.h */,
+ 6561B0F314867BCF0096B9CB /* statistics.c */,
+ 6561B0F614867BD70096B9CB /* statistics.h */,
65E7DFEF14600AD300095D43 /* qed.c */,
);
path = qed;
@@ -282,6 +287,7 @@
65E7E01414600F7B00095D43 /* qed.c in Sources */,
65AC1ABA14698B2B00336FFC /* leapfrog2.c in Sources */,
65AC1ABE14698B4400336FFC /* rec_lf_integrator.c in Sources */,
+ 6561B0F414867BCF0096B9CB /* statistics.c in Sources */,
);
runOnlyForDeploymentPostprocessing = 0;
};
View
70 qed/qed.c
@@ -5,6 +5,7 @@
#include <getopt.h>
#include <string.h>
+#include "statistics.h"
#include "lattice.h"
#include "hmc.h"
#include "fields.h"
@@ -137,11 +138,6 @@ int main(int argc, char **argv)
int total_updates = 0; //Total number of updates
/* Containers for the mean plaquette, Polyakov loop and chiral condensate */
/* and corresponding errors etc. */
- double mp = 0.0, dmp = 0.0, amp = 0.0, pl = 0.0, dpl = 0.0, apl = 0.0, cc = 0.0, dcc = 0.0, acc = 0.0;
- double emp, epl, ecc;
- // Containers for the mean value of Delta H and exp(-Delta H)
- double dh = 0, ddh = 0, adh = 0, edh = 0;
- double expdh = 0, dexpdh = 0, aexpdh = 0, eexpdh = 0;
double total_cgiterations1 = 0, total_cgiterations2 = 0;
/* Initialize the random number generator */
rlxd_init(2, 123456);
@@ -192,7 +188,13 @@ int main(int argc, char **argv)
gauge_force_sum = 0;
PF1_force_sum = 0;
PF2_force_sum = 0;
- for(i=0; i<g_measurements; i++)
+
+ statistics_data mp_statistics; reset_statistics_data(&mp_statistics);
+ statistics_data pl_statistics; reset_statistics_data(&pl_statistics);
+ statistics_data cc_statistics; reset_statistics_data(&cc_statistics);
+ statistics_data dh_statistics; reset_statistics_data(&dh_statistics);
+ statistics_data expdh_statistics; reset_statistics_data(&expdh_statistics);
+ for(i = 0; i < g_measurements; i ++)
{
g_cgiterations1 = 0;
g_cgiterations2 = 0;
@@ -203,27 +205,25 @@ int main(int argc, char **argv)
accepted_cur += update();
accepted += accepted_cur;
/* Measurements */
- mp = mean_plaquette();
- pl = polyakov_loop();
- cc = chiral_condensate();
- dh = ham_old - ham;
- amp += mp; //Mean plaquette
- dmp += mp*mp; // ... and its dispersion - for error estimation
- apl += pl; //Polyakov loop
- dpl += pl*pl;
- acc += cc; //Chiral condensate
- dcc += cc*cc;
- adh += dh;
- ddh += dh * dh;
- expdh = exp(dh);
- aexpdh += expdh;
- dexpdh += expdh * expdh;
+ double mp = mean_plaquette();
+ double pl = polyakov_loop();
+ double cc = chiral_condensate();
+
+ for (int t = 0; t < X2; t ++)
+ printf("%f\n", pion_correlation_function(t));
+
+ double dh = ham_old - ham;
+ add_statistics_entry(&mp_statistics, mp);
+ add_statistics_entry(&pl_statistics, pl);
+ add_statistics_entry(&cc_statistics, cc);
+ add_statistics_entry(&dh_statistics, dh);
+ add_statistics_entry(&expdh_statistics, exp(dh));
total_cgiterations1 += g_cgiterations1;
total_cgiterations2 += g_cgiterations2;
printf("\t Step %04i,\t mp = %2.4lf,\t pl = %2.4lf,\t cc = %2.4lf,\t dh = %2.4lf,\tcg1 = %d,\tcg2 = %d,\tacc = %d\n", i, mp, pl, cc, -dh, g_cgiterations1, g_cgiterations2, accepted_cur);
- };
+ }
/* Some output for diagnostics */
printf("\n\n Algorithm configuration:\n");
@@ -246,28 +246,12 @@ int main(int argc, char **argv)
printf("\t Runtime / seconds: %2.2lf\n", (double)(clock() - clock_start)/(double)CLOCKS_PER_SEC);
/* Measurement results */
- amp = amp/(double)g_measurements;
- apl = apl/(double)g_measurements;
- acc = acc/(double)g_measurements;
- adh = adh/(double)g_measurements;
- aexpdh = aexpdh/(double)g_measurements;
- dmp = dmp/(double)g_measurements - amp*amp;
- dpl = dpl/(double)g_measurements - apl*apl;
- dcc = dcc/(double)g_measurements - acc*acc;
- ddh = ddh/(double)g_measurements - adh*adh;
- dexpdh = dexpdh/(double)g_measurements - aexpdh*acc;
- emp = sqrt(dmp/(double)(g_measurements - 1)); //Errors of corresponding mean values
- epl = sqrt(dpl/(double)(g_measurements - 1));
- ecc = sqrt(dcc/(double)(g_measurements - 1));
- edh = sqrt(ddh/(double)(g_measurements - 1));
- eexpdh = sqrt(dexpdh/(double)(g_measurements - 1));
-
printf("\n\n Mean values:\n");
- printf("\t Plaquette: %2.6lf +/- %2.6lf (e = %02i%%)\n", amp, emp, (int)(100.0*emp/amp));
- printf("\t Polyakov loop: %2.6lf +/- %2.6lf (e = %02i%%)\n", apl, epl, (int)(100.0*epl/apl));
- printf("\t Chiral Condensate: %2.6lf +/- %2.6lf (e = %02i%%)\n", acc, ecc, (int)(100.0*ecc/acc));
- printf("\t -Delta H: %2.6lf +/- %2.6lf (e = %02i%%)\n", adh, edh, (int)(100.0*edh/adh));
- printf("\t exp(-Delta H): %2.6lf +/- %2.6lf (e = %02i%%)\n", aexpdh, eexpdh, (int)(100.0*eexpdh/aexpdh));
+ print_statistics_data(&mp_statistics, "Plaquette:");
+ print_statistics_data(&pl_statistics, "Polyakov loop:");
+ print_statistics_data(&cc_statistics, "Chiral Condensate:");
+ print_statistics_data(&dh_statistics, "-Delta H:");
+ print_statistics_data(&expdh_statistics, "exp(-Delta H):");
printf("\n\n");
free(left1);
View
45 qed/statistics.c
@@ -0,0 +1,45 @@
+//
+// statistics.c
+// qed
+//
+// Created by Daniel Alm on 30.11.11.
+// Copyright (c) 2011 __MyCompanyName__. All rights reserved.
+//
+
+#include <stdio.h>
+
+#include <string.h>
+#include <math.h>
+
+#include "statistics.h"
+
+void reset_statistics_data(statistics_data *data)
+{
+ memset(data, 0, sizeof(statistics_data));
+}
+
+void add_statistics_entry(statistics_data *data, double entry)
+{
+ data->N += 1;
+ data->sum += entry;
+ data->square_sum += entry * entry;
+}
+
+void calculate_statistics_data(statistics_data *data)
+{
+ int N = data->N;
+ data->mean = data->sum / N;
+ data->error = sqrt((data->square_sum / N - data->mean * data->mean) / (N - 1));
+}
+
+void calculate_statistics_array(statistics_data *data, int array_size)
+{
+ for (int i = 0; i < array_size; i ++)
+ calculate_statistics_data(&data[i]);
+}
+
+void print_statistics_data(statistics_data *data, const char *name)
+{
+ calculate_statistics_data(data);
+ printf("\t %-28s%2.6lf +/- %2.6lf (e = %02i%%)\n", name, data->mean, data->error, (int)(100.0*data->error/data->mean));
+}
View
25 qed/statistics.h
@@ -0,0 +1,25 @@
+//
+// statistics.h
+// qed
+//
+// Created by Daniel Alm on 30.11.11.
+// Copyright (c) 2011 __MyCompanyName__. All rights reserved.
+//
+
+#ifndef qed_statistics_h
+#define qed_statistics_h
+
+typedef struct statistics_data {
+ int N;
+ double sum, square_sum, mean, error;
+} statistics_data;
+
+void reset_statistics_data(statistics_data *data);
+
+void add_statistics_entry(statistics_data *data, double entry);
+void calculate_statistics_data(statistics_data *data);
+void calculate_statistics_array(statistics_data *data, int array_size);
+
+void print_statistics_data(statistics_data *data, const char *name);
+
+#endif
Please sign in to comment.
Something went wrong with that request. Please try again.