muller.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <config.h>

#include <stdlib.h>
#include <stdio.h>
#include <inttypes.h>
#include <string.h>

#include <jwsc/base/error.h>
#include <jwsc/base/option.h>
#include <jwsc/math/constants.h>
#include <jwsc/stat/dynamics.h>

#include "sampler.h"


#define  NUM_OPTS_NO_ARG      2
#define  NUM_OPTS_WITH_ARG    15

#define  DEFAULT_SEED         289375
#define  DEFAULT_ITERATIONS   1000
#define  DEFAULT_HYBRID_LEAP_ITER 100
#define  DEFAULT_START_X      0
#define  DEFAULT_START_Y      0
#define  DEFAULT_X_SIGMA      0.1
#define  DEFAULT_X_MIN       -3.0
#define  DEFAULT_X_MAX        1.5
#define  DEFAULT_Y_SIGMA      0.1
#define  DEFAULT_Y_MIN       -1.0
#define  DEFAULT_Y_MAX        3.0
#define  DEFAULT_DELTA_T      0.01
#define  DEFAULT_ALPHA        0.9
#define  DEFAULT_KICK         0
#define  DEFAULT_HYPER_GAMMA  1.0e-5
#define  DEFAULT_HYPER_ETA    1.0e-2
#define  DEFAULT_HYPER_BIAS   150.0
#define  DEFAULT_HYPER_SCALE  0.25
#define  DEFAULT_SAMPLER      MH


typedef enum
{
    MH,
    LANGE,
    HYPER,
    STOCH,
    HYBRID
}
Sampler_type;


static Option_no_arg opts_no_arg[NUM_OPTS_NO_ARG];

static Option_with_arg opts_with_arg[NUM_OPTS_WITH_ARG];
    
static uint32_t seed = DEFAULT_SEED;

static uint32_t iterations = DEFAULT_ITERATIONS;

static uint32_t hybrid_leap_iter = DEFAULT_HYBRID_LEAP_ITER;

static double start_x = DEFAULT_START_X;

static double start_y = DEFAULT_START_Y;

static double x_sigma = DEFAULT_X_SIGMA;

static double x_min = DEFAULT_X_MIN;

static double x_max = DEFAULT_X_MAX;

static double y_sigma = DEFAULT_Y_SIGMA;

static double y_min = DEFAULT_Y_MIN;

static double y_max = DEFAULT_Y_MAX;

static double delta_t = DEFAULT_DELTA_T;

static double alpha = DEFAULT_ALPHA;

static uint32_t kick = DEFAULT_KICK;

static double hyper_bias = DEFAULT_HYPER_BIAS;

static double hyper_scale = DEFAULT_HYPER_SCALE;

static double hyper_gamma = DEFAULT_HYPER_GAMMA;

static double hyper_eta = DEFAULT_HYPER_ETA;

static Sampler_type sampler = DEFAULT_SAMPLER;



static void print_usage(void)
{
    fprintf(stderr, "usage: muller OPTIONS\n");
    fprintf(stderr, "where OPTIONS is one or more of the following:\n");
    print_options(stderr, 25, NUM_OPTS_NO_ARG, opts_no_arg, NUM_OPTS_WITH_ARG,
            opts_with_arg);
}

static Error* process_help_opt(void)
{
    print_usage();
    exit(EXIT_SUCCESS);
}

static Error* process_version_opt(void)
{
    fprintf(stderr, "%s\n", MULLER_PACKAGE_STRING);
    exit(EXIT_SUCCESS);
}

static Error* process_seed_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'seed' requires an argument");
    }
    if (sscanf(arg, "%u", &seed) != 1)
    {
        return JWSC_EARG("Option 'seed' requires an argument");
    }
    return NULL;
}

static Error* process_iterations_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'iterations' requires an argument");
    }
    if (sscanf(arg, "%u", &iterations) != 1)
    {
        return JWSC_EARG("Option 'iterations' requires an argument");
    }
    return NULL;
}

static Error* process_hybrid_leap_iter_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'hybrid-leap-iter' requires an argument");
    }
    if (sscanf(arg, "%u", &hybrid_leap_iter) != 1)
    {
        return JWSC_EARG("Option 'hybrid-leap-iter' requires an argument");
    }
    return NULL;
}

static Error* process_start_x_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'start-x' requires an argument");
    }
    if (sscanf(arg, "%lf", &start_x) != 1)
    {
        return JWSC_EARG("Option 'start-x' must be > 0");
    }
    return NULL;
}

static Error* process_start_y_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'start-y' requires an argument");
    }
    if (sscanf(arg, "%lf", &start_y) != 1)
    {
        return JWSC_EARG("Option 'start-y' must be > 0");
    }
    return NULL;
}

static Error* process_delta_t_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'delta-t' requires an argument");
    }
    if (sscanf(arg, "%lf", &delta_t) != 1 || delta_t <= 0)
    {
        return JWSC_EARG("Option 'delta-t' must be > 0");
    }
    return NULL;
}

static Error* process_alpha_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'alpha' requires an argument");
    }
    if (sscanf(arg, "%lf", &alpha) != 1 || alpha < 0 || alpha > 1)
    {
        return JWSC_EARG("Option 'alpha' must be in [0,1]");
    }
    return NULL;
}

static Error* process_kick_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'kick' requires an argument");
    }
    if (sscanf(arg, "%u", &kick) != 1)
    {
        return JWSC_EARG("Option 'kick' requires an argument");
    }
    return NULL;
}

static Error* process_x_lim_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'x-lim' requires an argument");
    }
    if (sscanf(arg, "%lf,%lf", &x_min, &x_max) != 2)
    {
        return JWSC_EARG("Option 'x-lim' has format min,max");
    }
    if (x_max < x_min)
    {
        return JWSC_EARG("Option 'x-lim' max must be >= min");
    }
    return NULL;
}

static Error* process_x_sigma_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'x-sigma' requires an argument");
    }
    if (sscanf(arg, "%lf", &x_sigma) != 1 || x_sigma <= 1.0e-8)
    {
        return JWSC_EARG("Option 'x-sigma' must be > 0");
    }
    return NULL;
}

static Error* process_y_lim_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'y-lim' requires an argument");
    }
    if (sscanf(arg, "%lf,%lf", &y_min, &y_max) != 2)
    {
        return JWSC_EARG("Option 'y-lim' has format min,max");
    }
    if (y_max < y_min)
    {
        return JWSC_EARG("Option 'y-lim' max must be >= min");
    }
    return NULL;
}

static Error* process_y_sigma_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'y-sigma' requires an argument");
    }
    if (sscanf(arg, "%lf", &y_sigma) != 1 || y_sigma <= 1.0e-8)
    {
        return JWSC_EARG("Option 'y-sigma' must be > 0");
    }
    return NULL;
}

static Error* process_hyper_bias_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'hyper-bias' requires an argument");
    }
    if (sscanf(arg, "%lf", &hyper_bias) != 1)
    {
        return JWSC_EARG("Option 'hyper-bias' requires an argument");
    }
    return NULL;
}

static Error* process_hyper_scale_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'hyper-scale' requires an argument");
    }
    if (sscanf(arg, "%lf", &hyper_scale) != 1 || hyper_scale <= 1.0e-8)
    {
        return JWSC_EARG("Option 'hyper-scale' must be > 0");
    }
    return NULL;
}

static Error* process_sampler_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'sampler' requires an argument");
    }
    if (strncmp(arg, "mh", 2) == 0)
    {
        sampler = MH;
    }
    else if (strncmp(arg, "lange", 5) == 0)
    {
        sampler = LANGE;
    }
    else if (strncmp(arg, "hyper", 5) == 0)
    {
        sampler = HYPER;
    }
    else if (strncmp(arg, "stoch", 5) == 0)
    {
        sampler = STOCH;
    }
    else if (strncmp(arg, "hybrid", 6) == 0)
    {
        sampler = HYBRID;
    }
    else
    {
        return JWSC_EARG("Option 'sampler' must be one of {mh, lange, hyper, stoch, hybrid}");
    }
    return NULL;
}




int main(int argc, const char** argv)
{
    int      argi;
    uint32_t i;

    Params p;

    Error*   err;

    Option_string        l_name;
    Option_flag          s_name;
    Option_desc          desc;
    Option_func_no_arg   fnoarg;
    Option_func_with_arg farg;


    i      = 0;
    l_name = "help";
    s_name = 'h';
    desc   = "Prints program usage.";
    fnoarg = process_help_opt;
    init_option_no_arg(&(opts_no_arg[i++]), l_name, s_name, desc, fnoarg);

    l_name = "version";
    s_name = 'v';
    desc   = "Prints program version.";
    fnoarg = process_version_opt;
    init_option_no_arg(&(opts_no_arg[i++]), l_name, s_name, desc, fnoarg);

    i      = 0;
    l_name = "seed";
    s_name = 0;
    desc   = "Random seed.";
    farg   = process_seed_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "iterations";
    s_name = 0;
    desc   = "Number of sampling iterations.";
    farg   = process_iterations_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "hybrid-leap-iter";
    s_name = 0;
    desc   = "Number of hybrid Monte Carlo leapfrog iterations.";
    farg   = process_hybrid_leap_iter_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "start-x";
    s_name = 0;
    desc   = "Starting x-coord position.";
    farg   = process_start_x_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "start-y";
    s_name = 0;
    desc   = "Starting y-coord position.";
    farg   = process_start_y_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "delta-t";
    s_name = 0;
    desc   = "Stochastic dynamics integration step-size for Langevin Monte Carlo.";
    farg   = process_delta_t_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "alpha";
    s_name = 0;
    desc   = "Stochastic dynamics step-size for Langevin Monte Carlo.";
    farg   = process_alpha_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "kick";
    s_name = 0;
    desc   = "Stochastic dynamics frequency to kick the particle by resetting the momenta. Zero disables any kicking. Must be > 0.";
    farg   = process_kick_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "x-sigma";
    s_name = 0;
    desc   = "Gaussian sigma for MH proposal of x. Must be > 0.";
    farg   = process_x_sigma_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "x-lim";
    s_name = 0;
    desc   = "Limits over the range of x. Has format min,max";
    farg   = process_x_lim_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "y-sigma";
    s_name = 0;
    desc   = "Gaussian sigma for MH proposal of y. Must be > 0.";
    farg   = process_y_sigma_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "y-lim";
    s_name = 0;
    desc   = "Limits over the range of y. Has format min,max.";
    farg   = process_y_lim_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "sampler";
    s_name = 0;
    desc   = "Type of sampler to use. Must be one of {mh, lange, hyper, stoch, hybrid}.";
    farg   = process_sampler_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "hyper-bias";
    s_name = 0;
    desc   = "Bias parameter for hyper sampling.";
    farg   = process_hyper_bias_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "hyper-scale";
    s_name = 0;
    desc   = "Scale parameter for hyper sampling.";
    farg   = process_hyper_scale_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);


    if ((err = process_options(argc, argv, &argi, NUM_OPTS_NO_ARG, opts_no_arg,
                    NUM_OPTS_WITH_ARG, opts_with_arg)) != NULL)
    {
        print_error_msg_exit("muller", err->msg);
    }

    srand(seed); rand(); rand();

    if (start_x == DEFAULT_START_X)
    {
        start_x = sample_uniform_pdf_d(x_min, x_max);
    }
    if (start_y == DEFAULT_START_Y)
    {
        start_y = sample_uniform_pdf_d(y_min, y_max);
    }

    p.x = start_x;
    p.y = start_y;

    if (sampler == MH)
    {
        if ((err = metropolis_hastings(iterations, start_x, start_y, x_sigma,
                        x_min, x_max, y_sigma, y_min, y_max)) != NULL)
        {
            print_error_msg_exit("muller", err->msg);
        }
    }
    else if (sampler == LANGE)
    {
        if ((err = langevin_1_d(iterations, delta_t, &p, get_pt_from_params,
                        grad_energy_func, langevin_accept_sample)) 
                != NULL)
        {
            print_error_msg_exit("muller", err->msg);
        }
    }
    else if (sampler == HYPER)
    {
        if ((err = langevin_hyperdynamics_1_d(iterations, delta_t, hyper_gamma,
                        hyper_eta, hyper_bias, hyper_scale, &p,
                        get_pt_from_params, set_params_from_pt, energy_func,
                        grad_energy_func, hyper_accept_sample)) != NULL)
        {
            print_error_msg_exit("muller", err->msg);
        }
    }
    else if (sampler == STOCH)
    {
        if ((err = stochastic_dynamics_1_d(iterations, delta_t, alpha, kick, &p,
                        get_pt_from_params, grad_energy_func,
                        stochastic_dynamics_accept_sample)) != NULL)
        {
            print_error_msg_exit("muller", err->msg);
        }
    }
    else if (sampler == HYBRID)
    {
        if ((err = hybrid_monte_carlo_1_d(iterations, hybrid_leap_iter,
                        delta_t, alpha, kick, &p, get_pt_from_params,
                        set_params_from_pt, energy_func, grad_energy_func,
                        hybrid_mc_accept_sample, NULL)) != NULL)
        {
            print_error_msg_exit("muller", err->msg);
        }
    }

    if (get_num_unhandled_errors() > 0)
    {
        print_error_msg_exit("muller", "Unhandled errors");
    }

    return EXIT_SUCCESS;
}