From 9303f4f3793ab5f158cee63e50a7c172de69989a Mon Sep 17 00:00:00 2001 From: rnhmjoj Date: Mon, 18 May 2020 00:32:32 +0200 Subject: [PATCH] ex-6 compute EMD of the convolution --- ex-6/dist-plot.py | 22 +++++++++++----- ex-6/test.c | 65 +++++++++++++++++++++++++---------------------- 2 files changed, 50 insertions(+), 37 deletions(-) diff --git a/ex-6/dist-plot.py b/ex-6/dist-plot.py index 56bedff..99e707a 100755 --- a/ex-6/dist-plot.py +++ b/ex-6/dist-plot.py @@ -8,21 +8,29 @@ rcParams['font.size'] = 12 n = int(input()) a, b, f = loadtxt(sys.stdin, unpack=True) -subplot(121) +subplot(131) title('FFT') -xlabel('EDM distance') +xlabel('EDM') ylabel('counts') hist(a[:n], insert(b[:n], 0, a[0]), weights=f[:n], histtype='stepfilled', color='#e3c5ca', edgecolor='#92182b') -ticklabel_format(style='sci', axis='x', scilimits=(0, 0)) +ticklabel_format(style='sci', axis='x', scilimits=(0, 0), useMathText=True) -subplot(122) +subplot(132) title('RL') -xlabel('EDM distance') +xlabel('EDM') ylabel('counts') -hist(a[n:], insert(b[n:], 0, a[n]), weights=f[n:], +hist(a[n:2*n], insert(b[n:2*n], 0, a[n]), weights=f[n:2*n], histtype='stepfilled', color='#e3c5ca', edgecolor='#92182b') -ticklabel_format(style='sci', axis='x', scilimits=(0, 0)) +ticklabel_format(style='sci', axis='x', scilimits=(0, 0), useMathText=True) + +subplot(133) +title('convolution') +xlabel('EDM') +ylabel('counts') +hist(a[2*n:], insert(b[2*n:], 0, a[2*n]), weights=f[2*n:], + histtype='stepfilled', color='#e3c5ca', edgecolor='#92182b') +ticklabel_format(style='sci', axis='x', scilimits=(0, 0), useMathText=True) tight_layout() show() diff --git a/ex-6/test.c b/ex-6/test.c index 5824c95..85eeb93 100644 --- a/ex-6/test.c +++ b/ex-6/test.c @@ -21,14 +21,6 @@ struct options { double noise; }; - -/* A pair of `double`s */ -typedef struct { - double a; - double b; -} pair_t; - - int show_help(char **argv) { fprintf(stderr, "Usage: %s -[hexbsrn]\n", argv[0]); fprintf(stderr, " -h\t\tShow this message.\n"); @@ -61,7 +53,7 @@ int show_help(char **argv) { * The function returns a pair of the FFT and RL * distances, in this order. */ -pair_t experiment( +double *experiment( struct options opts, gsl_rng *r, gsl_histogram *kernel) { @@ -103,9 +95,10 @@ pair_t experiment( * histogram and store add each one to the respective * distance histogram. */ - pair_t dist; - dist.a = emd_between(hist, fft_clean); - dist.b = emd_between(hist, rl_clean); + static double dist[3]; + dist[0] = emd_between(hist, fft_clean); + dist[1] = emd_between(hist, rl_clean); + dist[2] = emd_between(hist, conv); // free memory gsl_histogram_free(hist); @@ -149,13 +142,15 @@ int main(int argc, char **argv) { /* Performs experiments and record the deconvolution * EMD distances to the original sample. */ - double* fft_dist = calloc(opts.num_exps, sizeof(double)); - double* rl_dist = calloc(opts.num_exps, sizeof(double)); + double* fft_dist = calloc(opts.num_exps, sizeof(double)); + double* rl_dist = calloc(opts.num_exps, sizeof(double)); + double* conv_dist = calloc(opts.num_exps, sizeof(double)); for (size_t i = 0; i < opts.num_exps; i++) { - pair_t dist = experiment(opts, r, kernel); - fft_dist[i] = dist.a; - rl_dist[i] = dist.b; + double *dist = experiment(opts, r, kernel); + fft_dist[i] = dist[0]; + rl_dist[i] = dist[1]; + conv_dist[i] = dist[2]; } /* Compute some statistics of the EMD */ @@ -175,18 +170,21 @@ int main(int argc, char **argv) { double conv_skew = gsl_stats_skew_m_sd(conv_dist, 1, opts.num_exps, conv_mean, conv_stdev); double conv_min, conv_max; gsl_stats_minmax(&conv_min, &conv_max, conv_dist, 1, opts.num_exps); - /* Create EMD distance histograms. + /* Create EDM distance histograms. * Since the distance depends wildly on the noise we can't * set a fixed range and therefore use the above values. */ - gsl_histogram *fft_hist = gsl_histogram_alloc(sqrt(opts.num_exps)); - gsl_histogram *rl_hist = gsl_histogram_alloc(sqrt(opts.num_exps)); - gsl_histogram_set_ranges_uniform(fft_hist, fft_min, fft_max); - gsl_histogram_set_ranges_uniform(rl_hist, rl_min, rl_max); + gsl_histogram *fft_hist = gsl_histogram_alloc(sqrt(opts.num_exps)); + gsl_histogram *rl_hist = gsl_histogram_alloc(sqrt(opts.num_exps)); + gsl_histogram *conv_hist = gsl_histogram_alloc(sqrt(opts.num_exps)); + gsl_histogram_set_ranges_uniform( fft_hist, fft_min, fft_max); + gsl_histogram_set_ranges_uniform( rl_hist, rl_min, rl_max); + gsl_histogram_set_ranges_uniform(conv_hist, conv_min, conv_max); for (size_t i = 0; i < opts.num_exps; i++) { - gsl_histogram_increment(fft_hist, fft_dist[i]); - gsl_histogram_increment( rl_hist, rl_dist[i]); + gsl_histogram_increment( fft_hist, fft_dist[i]); + gsl_histogram_increment( rl_hist, rl_dist[i]); + gsl_histogram_increment(conv_hist, conv_dist[i]); } // print results to stderr @@ -199,20 +197,27 @@ int main(int argc, char **argv) { fprintf(stderr, "- mean: %.2e\n" "- stdev: %.1e\n" "- skew: %.2f\n", rl_mean, rl_stdev, rl_skew); + fprintf(stderr, "\n## convolution\n"); + fprintf(stderr, "- mean: %.2e\n" + "- stdev: %.1e\n" + "- skew: %.2f\n", conv_mean, conv_stdev, conv_skew); // print histograms to stdout fprintf(stderr, "\n# EMD histogram\n"); fprintf(stdout, "%d\n", (int)sqrt(opts.num_exps)); - gsl_histogram_fprintf(stdout, fft_hist, "%g", "%g"); - gsl_histogram_fprintf(stdout, rl_hist, "%g", "%g"); + gsl_histogram_fprintf(stdout, fft_hist, "%.10g", "%.10g"); + gsl_histogram_fprintf(stdout, rl_hist, "%.10g", "%.10g"); + gsl_histogram_fprintf(stdout, conv_hist, "%.10g", "%.10g"); // free memory gsl_rng_free(r); gsl_histogram_free(kernel); - gsl_histogram_free(fft_hist); - gsl_histogram_free(rl_hist); - free(fft_dist); - free(rl_dist); + gsl_histogram_free( fft_hist); + gsl_histogram_free( rl_hist); + gsl_histogram_free(conv_hist); + free( fft_dist); + free( rl_dist); + free(conv_dist); return EXIT_SUCCESS; }