gaussian_bench.c

Go to the documentation of this file.

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "basic.h"
#include "strategy_gaussian_conv.h"
#include "filter_util.h"

#define OUTPUT_FILE     "impulse.txt"

#ifndef M_SQRT2PI

#define M_SQRT2PI   2.50662827463100050241576528481104525
#endif


void print_usage()
{
    puts("Gaussian benchmark, P. Getreuer 2012-2013");
#ifdef NUM_SINGLE
    puts("Configuration: single-precision computation");
#else
    puts("Configuration: double-precision computation");
#endif
    puts("\nSyntax: gaussian_bench [bench type] [options] [output]\n");
    puts("Bench type:");
    puts("   speed         measure computation time");
    puts("   accuracy      measure L^infty operator norm error");
    puts("   impulse       impulse response, written to " OUTPUT_FILE "\n");
    puts("Options:");
    puts("   -a <algo>     algorithm to use, choices are");
    puts("                 fir     FIR approximation, tol = kernel accuracy");
    puts("                 dct     DCT-based convolution");
    puts("                 box     box filtering, K = # passes");
    puts("                 sii     stacked integral images, K = # boxes");
    puts("                 am_orig Alvarez-Mazorra original method,");
    puts("                         K = # passes, tol = boundary accuracy");
    puts("                 am      Alvarez-Mazorra using regression on q,");
    puts("                         K = # passes, tol = boundary accuracy");
    puts("                 deriche Deriche recursive filtering,");
    puts("                         K = order, tol = boundary accuracy");
    puts("                 vyv     Vliet-Young-Verbeek recursive filtering,");
    puts("                         K = order, tol = boundary accuracy");
    puts("   -s <number>   sigma, standard deviation of the Gaussian");
    puts("   -K <number>   specifies number of steps (box, sii, am)");
    puts("   -t <number>   accuracy tolerance (fir, am, deriche, vyv)");
    puts("   -N <number>   signal length\n");
    puts("   -r <number>   (speed bench) number of runs");
    puts("   -n <number>   (impulse bench) position of the impulse\n");
}

typedef struct
{
    const char *bench_type;     
    long N;                     
    long num_runs;              
    long n0;                    
    const char *algo;           
    double sigma;               
    int K;                      
    double tol;                 
} program_params;

int parse_params(program_params *param, int argc, char **argv);

int speed_test(program_params p, num *output, num *input)
{
    gconv *g = NULL;
    unsigned long time_start, time_stop;
    long run;
    
    if (!(g = gconv_plan(output, input, p.N, 1,
        p.algo, p.sigma, p.K, p.tol)))
        return 0;
    
    time_start = millisecond_timer();
    
    for (run = 0; run < p.num_runs; ++run)
        gconv_execute(g);
        
    time_stop = millisecond_timer();
    printf("%.5e\n",
        ((double)(time_stop - time_start)) / p.num_runs);
    gconv_free(g);
    return 1;
}

void make_impulse_signal(num *signal, long N, long n0)
{
    long n;
    
    for (n = 0; n < N; ++n)
        signal[n] = 0;
    
    signal[n0] = 1;
}

int accuracy_test(program_params p, num *output, num *input)
{
    double *error_sums = NULL;
    num *output0 = NULL;
    gconv *g0 = NULL, *g = NULL;
    double linf_norm = 0.0;
    long m, n;
    int success = 0;
        
    if (!(error_sums = (double *)malloc(sizeof(double) * p.N))
        || !(output0 = (num *)malloc(sizeof(num) * p.N))
        || !(g0 = gconv_plan(output0, input, p.N, 1,
            "fir", p.sigma, p.K, 1e-15))
        || !(g = gconv_plan(output, input, p.N, 1,
            p.algo, p.sigma, p.K, p.tol)))
        goto fail;
    
    for (n = 0; n < p.N; ++n)
        error_sums[n] = 0.0;
    
    for (n = 0; n < p.N; ++n)
    {
        make_impulse_signal(input, p.N, n);
        gconv_execute(g0);
        gconv_execute(g);
        
        for (m = 0; m < p.N; ++m)
            error_sums[m] += fabs((double)output0[m] - (double)output[m]);
    }
    
    for (n = 0; n < p.N; ++n)
        if (error_sums[n] > linf_norm)
            linf_norm = error_sums[n];
    
    printf("%.8e\n", linf_norm);
    success = 1;
fail:
    gconv_free(g0);
    gconv_free(g);
    if (output0)
        free(output0);
    if (error_sums)
        free(error_sums);
    return success;
}

int write_output(const char *filename, num *output, num *exact, long N)
{
    FILE *f;
    long n;
    
    if (!(f = fopen(filename, "wt")))
        return 0;
    
    fprintf(f, "# n\toutput value\texact value\n");
    
    for (n = 0; n < N; ++n)
        fprintf(f, "%ld\t%.16e\t%.16e\n", n, output[n], exact[n]);
    
    fclose(f);
    return 1;
}

int impulse_test(program_params p, num *output, num *input)
{
    num *exact = NULL;
    gconv *g = NULL, *g_exact = NULL;
    int success = 0;
    
    if (!(exact = (num *)malloc(sizeof(num) * p.N))
        || !(g = gconv_plan(output, input, p.N, 1,
            p.algo, p.sigma, p.K, p.tol))
        || !(g_exact = gconv_plan(exact, input, p.N, 1,
            "fir", p.sigma, p.K, 1e-15)))
        goto fail;
    
    make_impulse_signal(input, p.N, p.n0);
    gconv_execute(g);
    gconv_execute(g_exact);
    
    if (!write_output(OUTPUT_FILE, output, exact, p.N))
        goto fail;
    
    success = 1;
fail:
    gconv_free(g_exact);
    gconv_free(g);
    if (exact)
        free(exact);
    return success;
}

int main(int argc, char **argv)
{
    program_params param;
    num *input = NULL;
    num *output = NULL;
    int success = 0;
    
    if (!parse_params(&param, argc, argv))
        return 1;
    
    /* Create input signal */
    if (!(input = (num *)malloc(sizeof(num) * param.N))
        || !(output = (num *)malloc(sizeof(num) * param.N)))
    {
        fprintf(stderr, "Out of memory\n");
        goto fail;
    }
    
    if (!strcmp(param.bench_type, "speed"))
    {
        if (!speed_test(param, output, input))
            goto fail;
    }
    else if (!strcmp(param.bench_type, "accuracy"))
    {
        if (!accuracy_test(param, output, input))
            goto fail;
    }
    else if (!strcmp(param.bench_type, "impulse"))
    {
        if (!impulse_test(param, output, input))
            goto fail;
    }
    else
    {
        fprintf(stderr, "Invalid bench type \"%s\"\n", param.bench_type);
        goto fail;
    }
    
    success = 1;
fail:
    if (!success)
        fprintf(stderr, "Bench %s failed\n", param.bench_type);

    if (output)
        free(output);
    if (input)
        free(input);
    return !success;
}

int parse_params(program_params *param, int argc, char **argv)
{
    static const char *default_algo = (const char *)"exact";
    char *option_string;
    char option_char;
    int i;

    if (argc < 2)
    {
        print_usage();
        return 0;
    }

    param->bench_type = argv[1];
    
    /* Set parameter defaults. */
    param->N = 100;
    param->n0 = -1;
    param->num_runs = 10;
    param->sigma = 5;
    param->algo = default_algo;
    param->K = 3;
    param->tol = 1e-2;
    
    for (i = 2; i < argc;)
    {
        if (argv[i] && argv[i][0] == '-')
        {
            if ((option_char = argv[i][1]) == 0)
            {
                fprintf(stderr, "Invalid parameter format.\n");
                return 0;
            }

            if (argv[i][2])
                option_string = &argv[i][2];
            else if (++i < argc)
                option_string = argv[i];
            else
            {
                fprintf(stderr, "Invalid parameter format.\n");
                return 0;
            }
            
            switch (option_char)
            {
            case 'a':   /* Read algorithm. */
                param->algo = option_string;
                break;
            case 's':   /* Read sigma parameter. */
                param->sigma = atof(option_string);
                
                if (param->sigma < 0)
                {
                    fprintf(stderr, "sigma must be positive.\n");
                    return 0;
                }
                break;
            case 'K':   /* Read number of steps. */
                param->K = atoi(option_string);
                
                if (param->K <= 0)
                {
                    fprintf(stderr, "K must be positive.\n");
                    return 0;
                }
                break;
            case 't':   /* Read tolerance. */
                param->tol = atof(option_string);
                
                if (param->tol < 0)
                {
                    fprintf(stderr, "Tolerance must be positive.\n");
                    return 0;
                }
                break;
            case 'N':   /* Input length. */
                param->N = atoi(option_string);
                
                if (param->N < 0)
                {
                    fprintf(stderr, "Signal length must be positive.\n");
                    return 0;
                }
                break;
            case 'r':   /* Read number of runs. */
                param->num_runs = (long)atof(option_string);
                
                if (param->num_runs <= 0)
                {
                    fprintf(stderr, "Number of runs must be positive.\n");
                    return 0;
                }
                break;
            case 'n':   /* Impulse position. */
                param->n0 = atoi(option_string);
                break;
            case '-':
                print_usage();
                return 0;
            default:
                if (isprint(option_char))
                    fprintf(stderr, "Unknown option \"-%c\".\n", option_char);
                else
                    fprintf(stderr, "Unknown option.\n");
            
                return 0;
            }
            
            i++;
        }
        else
            i++;
    }
    
    if (param->n0 < 0 || param->n0 >= param->N)
        param->n0 = (param->N + 1) / 2;
    
    return 1;
}