diff --git a/ex-5/lib.c b/ex-5/lib.c new file mode 100644 index 0000000..92528c4 --- /dev/null +++ b/ex-5/lib.c @@ -0,0 +1,86 @@ +#include +#include +#include +#include "lib.h" + + +// Wrapper for the gsl_function structure +// +double function (double * x, size_t dim, void * params) + { + return exp(x[0]); + } + +/////////////////////////////////////////////////////////////////////////////// + +// Results printer. +// +void results (size_t calls, double result, double error, double chi) + { + if (calls != 0) + { + printf ("%6.0e | ", (double)calls); + } + printf ("%5f | ", result); + printf ("%5f | ", error); + if (chi != 0) + { + printf ("%5f", chi); + } + } + +/////////////////////////////////////////////////////////////////////////////// + +// Perform a fit in order to compare the data with the expected function: +// +// y = a⋅x^b → ln(y) = ln(a) + b⋅ln(x) → Y = A + B⋅X +// +// with: +// +// A = ln(a) → a = e^A +// B = b → b = B +// +// using a linear regression. For b, the results is hence compared with the +// expected one: +// +// b_exp = - 0.5 +// +void fit (struct bag full_bag, double* p, + double* a, double* a_err, + double* b, double* b_err) + { + // Expected value. + // + double b_exp = -0.5; + + // Parse arguments. + // + double size = full_bag.size; + double* x = full_bag.pokets.x; + double* y = full_bag.pokets.y; + for (size_t i = 0; i < size; i++) + { + x[i] = log(x[i]); + y[i] = log(y[i]); + } + + // Do fit. + // + double A, B, A_err, B_err, AB_cov, sum2; + gsl_fit_linear(x, 1, y, 1, size, &A, &B, &A_err, &AB_cov, &B_err, &sum2); + + // Parse results. + // + A_err = sqrt(A_err); + B_err = sqrt(B_err); + *a = exp(A); + *b = B; + *a_err = *a * A_err; + *b_err = B_err; + + // Check compatibility with expected values. + // + double t = fabs(*b - b_exp)/ *b_err; + *p = 1 - erf(t/sqrt(2)); + } + diff --git a/ex-5/lib.h b/ex-5/lib.h new file mode 100644 index 0000000..6c7bf8f --- /dev/null +++ b/ex-5/lib.h @@ -0,0 +1,35 @@ +#pragma once + + +// Subtruct for the results storage. +// +struct poket + { + double * x; // Vector with the numbers of function calls. + double * y; // Vector with Plain MC error estimates. + }; + +// Struct for the results storage. +// +struct bag + { + struct poket pokets; // Values. + size_t size; // Quantity of values. + }; + +/////////////////////////////////////////////////////////////////////////////// + +// Wrapper for the gsl_function structure. +// +double function (double * x, size_t dim, void * params); + +// Results printer. +// +void results (size_t calls, double result, double error, double chi); + +// Do a linear minimization fitting the model y = a⋅x^b. The p-value of the +// compatibility test of b with the expected value -0.5 is returned in p. +// +void fit (struct bag full_bag, double* p, + double* a, double* a_err, + double* b, double* b_err); diff --git a/ex-5/main.c b/ex-5/main.c index 115cb1c..4b58526 100644 --- a/ex-5/main.c +++ b/ex-5/main.c @@ -2,43 +2,34 @@ #include #include #include +#include #include - - -// exp() wrapper for the gsl_function structure -// -double function (double * x, size_t dim, void * params) - { - return exp(x[0]); - } - - -// Function which prints out the results when called. -// -void results (size_t calls, double result, double error, double chi) - { - if (calls != 0) - { - printf ("%6.0e | ", (double)calls); - } - printf ("%5f | ", result); - printf ("%5f | ", error); - if (chi != 0) - { - printf ("%5f", chi); - } - } +#include "lib.h" int main(int argc, char** argv) { // Some useful variables. // - size_t dims = 1; // Integral dimension - double lower[1] = {0}; // Integration lower limit - double upper[1] = {1}; // Integration upper limit - size_t calls = 50; // Initial number of function calls - double integral, error; + size_t dims = 1; // Integral dimension + double lower[1] = {0}; // Integration range lower limit + double upper[1] = {1}; // Integration range upper limit + double integral, error; // Result and error + + // An integral is estimated with a number c_i of function calls. + // Different number of function calls are tested. + // + size_t c_0 = 50; // Initial number of function calls + size_t c_f = 50000000; // Final number of function calls + size_t factor = 10; // c_(i+1) = c_i * factor + size_t size = round(1 + log(c_f/c_0)/log(factor)); // Number of integrations + + // A fit will be performed. This struct is needed to accomplish it. + // + struct bag full_bag; + full_bag.pokets.x = calloc(size, sizeof(double)); + full_bag.pokets.y = calloc(size, sizeof(double)); + full_bag.size = size; // Initialize an RNG. // @@ -52,23 +43,18 @@ int main(int argc, char** argv) expo.dim = 1; expo.params = NULL; - // Print the results table header + // Print the results table header. + // printf(" Calls | Plain | Error | MISER |" " Error | VEGAS | Error | χ²\n"); printf(" ------|----------|----------|----------|" "----------|----------|----------|---------\n"); - // Compute the integral for the following number of function calls: - // 50 - // 500 - // 5'000 - // 50'000 - // 500'000 - // 5'000'000 - // 50'000'000 - // with the three MC methods: plain MC, MISER and VEGAS. + // Compute the integral with the three MC methods: plain MC, MISER and VEGAS. // - while (calls <= 5000000) + size_t calls = c_0; + size_t i = 0; + while (calls <= c_f) { // Plain MC gsl_monte_plain_state *sMC = gsl_monte_plain_alloc (dims); @@ -77,6 +63,11 @@ int main(int argc, char** argv) gsl_monte_plain_free(sMC); results(calls, integral, error, 0); + // Update the struct for the fit. + // + full_bag.pokets.x[i] = calls; + full_bag.pokets.y[i++] = error; + // MISER gsl_monte_miser_state *sMI = gsl_monte_miser_alloc (dims); gsl_monte_miser_integrate (&expo, lower, upper, dims, calls, @@ -106,13 +97,28 @@ int main(int argc, char** argv) double chi = sVE->chisq; results(0, integral, error, chi); - // Update function calls + // printf ("\n"); - calls = calls*10; + calls = calls*factor; } + // Do a fit of the Plain MC errors. + // + double p, a, a_err, b, b_err; + fit(full_bag, &p, &a, &a_err, &b, &b_err); + + // Print the fit results. + // + fprintf (stderr, "\n## Fit results:\n\n"); + fprintf (stderr, "a = %.5f\t", a); + fprintf (stderr, "δa = %.5f\n", a_err); + fprintf (stderr, "b = %.5f\t", b); + fprintf (stderr, "δb = %.5f\n", b_err); + fprintf (stderr, "p-value = %.5f\n", p); + // Free memory. + // gsl_rng_free(r); return EXIT_SUCCESS; diff --git a/ex-5/plot.py b/ex-5/plot.py index d581431..7ca0ba4 100755 --- a/ex-5/plot.py +++ b/ex-5/plot.py @@ -4,26 +4,29 @@ import matplotlib.pyplot as plt import numpy as np import sys + table = np.loadtxt(sys.stdin, unpack=True, skiprows=2, delimiter='|') calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE, chi = table exact = 1.7182818285 -plt.rcParams['font.size'] = 17 - plt.figure() +# plt.figure(figsize=(5, 3)) +# plt.rcParams['font.size'] = 8 + plt.title('Plain MC', loc='right') plt.ylabel('$I^{oss}$') plt.xlabel('calls') plt.xscale('log') -plt.axhline(y=exact, color='#c556ea', linestyle='-', - label='Exact value') -plt.errorbar(calls, I_MC, linestyle='', marker='o', - yerr=σ_MC, color='#92182b', label='Plain MC') -plt.errorbar(calls, I_MI, linestyle='', marker='o', - yerr=σ_MI, color='black', label='MISER') -plt.errorbar(calls, I_VE, linestyle='', marker='o', - yerr=σ_VE, color='gray', label='VEGAS') - +plt.axhline(y=exact, color='#c556ea', linestyle='-', label='Exact value') +plt.errorbar(calls, I_MC, linestyle='', marker='o', yerr=σ_MC, + color='#92182b', label='Plain MC') +plt.errorbar(calls, I_MI, linestyle='', marker='o', yerr=σ_MI, + color='black', label='MISER') +plt.errorbar(calls, I_VE, linestyle='', marker='o', yerr=σ_VE, + color='gray', label='VEGAS') plt.legend() + +# plt.tight_layout() +# plt.savefig('notes/images/5-test.pdf') plt.show() diff --git a/ex-5/plots/fit_plot.py b/ex-5/plots/fit_plot.py new file mode 100755 index 0000000..a2d40b0 --- /dev/null +++ b/ex-5/plots/fit_plot.py @@ -0,0 +1,31 @@ +#!/usr/bin/env python + +import matplotlib.pyplot as plt +import numpy as np +import sys + + +table = np.loadtxt(sys.stdin, unpack=True, skiprows=2, delimiter='|') +calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE, chi = table +exact = 1.7182818285 + +plt.figure() +# plt.figure(figsize=(5, 3)) +# plt.rcParams['font.size'] = 8 + +plt.title('Integral uncertainty', loc='right') +plt.ylabel('$\sigma_I$') +plt.xlabel('calls') +plt.xscale('log') +plt.yscale('log') + +plt.plot(calls, σ_MC, linestyle='', marker='o', + color='#92182b', label='$\sigma_I$') +x = np.logspace(1, 8, 5) +a = 0.48734 +b = -0.49938 +plt.plot(x, a*(x**b), color='gray') + +plt.tight_layout() +# plt.savefig('notes/images/5-fit.pdf') +plt.show() diff --git a/makefile b/makefile index 783009b..b0ff4d0 100644 --- a/makefile +++ b/makefile @@ -27,7 +27,7 @@ ex-4/bin/main: ex-4/main.c ex-4/lib.c $(CCOMPILE) ex-5: ex-5/bin/main -ex-5/bin/%: ex-5/%.c +ex-5/bin/%: ex-5/main.c ex-5/lib.c $(CCOMPILE) ex-6: ex-6/bin/main ex-6/bin/test diff --git a/notes/images/5-fit.pdf b/notes/images/5-fit.pdf new file mode 100644 index 0000000..d8bac2e Binary files /dev/null and b/notes/images/5-fit.pdf differ