sampler.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 <math.h>
#include <assert.h>
#include <inttypes.h>

#include <jwsc/base/error.h>
#include <jwsc/base/limits.h>
#include <jwsc/prob/pdf.h>
#include <jwsc/vector/vector.h>
#include <jwsc/vector/vector_math.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matrix/matrix_math.h>

#include "potential.h"
#include "sampler.h"


Error* metropolis_hastings
(
    uint32_t iterations,
    double   start_x,
    double   start_y,
    double   x_sigma,
    double   x_min,
    double   x_max,
    double   y_sigma,
    double   y_min,
    double   y_max
)
{
    uint32_t i;
    uint8_t  accepted;

    double log_alpha;
    double log_u;

    double x_accepted;
    double y_accepted;
    double V_accepted;

    double x_proposed;
    double y_proposed;
    double V_proposed;

    double q_accepted_lp;
    double q_proposed_lp;

    Error* err;

    x_accepted = start_x;
    y_accepted = start_y;
    x_proposed = x_accepted;
    y_proposed = y_accepted;

    if ((err = mullers_potential_d(&V_accepted, x_accepted, y_accepted)))
    {
        return err;
    }
    V_proposed = V_accepted;

    accepted = 1;

    for (i = 0; i < iterations; i++)
    {
        fprintf(stdout, "%7.4f  %7.4f  %8e  %7.4f  %7.4f  %8e  %1u\n",
                x_accepted,
                y_accepted,
                V_accepted, 
                x_proposed,
                y_proposed,
                V_proposed,
                accepted);

        x_proposed = sample_gaussian_pdf_d(x_accepted, x_sigma, x_min, x_max);
        y_proposed = sample_gaussian_pdf_d(y_accepted, y_sigma, y_min, y_max);

        if ((err = mullers_potential_d(&V_proposed, x_proposed, y_proposed)))
        {
            return err;
        }

        q_proposed_lp  = 0;
        q_proposed_lp += gaussian_pdf_d(x_accepted, x_sigma, x_proposed);
        q_proposed_lp += gaussian_pdf_d(y_accepted, y_sigma, y_proposed);

        q_accepted_lp  = 0;
        q_accepted_lp += gaussian_pdf_d(x_proposed, x_sigma, x_accepted);
        q_accepted_lp += gaussian_pdf_d(y_proposed, y_sigma, y_accepted);

        log_alpha  = 0;
        log_alpha += log((-V_proposed < JWSC_MIN_LOG_ARG) ? 
                         JWSC_MIN_LOG_ARG : -V_proposed);
        log_alpha -= log((-V_accepted < JWSC_MIN_LOG_ARG) ?
                         JWSC_MIN_LOG_ARG : -V_accepted);
        log_alpha += q_accepted_lp;
        log_alpha -= q_proposed_lp;

        log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));

        accepted = 0;

        if (log_u <= log_alpha)
        {
            accepted = 1;
            x_accepted = x_proposed;
            y_accepted = y_proposed;
            V_accepted = V_proposed;
        }
    }

    return NULL;
}


void get_pt_from_params(Vector_d** p_out, const void* params)
{
    create_vector_d(p_out, 2);
    (*p_out)->elts[0] = ((Params*)params)->x;
    (*p_out)->elts[1] = ((Params*)params)->y;
}

Error* set_params_from_pt(const Vector_d* p, void* params)
{
    ((Params*)params)->x = p->elts[0];
    ((Params*)params)->y = p->elts[1];

    return NULL;
}

Error* energy_func(double* energy_out, void* params)
{
    Error* err;

    if ((err = mullers_potential_d(energy_out, ((Params*)params)->x,
                    ((Params*)params)->y)))
    {
        return err;
    }

    return NULL;
}

Error* grad_energy_func(Vector_d** grad_out, void* params)
{
    Error* err;

    create_vector_d(grad_out, 2);

    if ((err = grad_mullers_potential_d(&((*grad_out)->elts[0]), 
                    &((*grad_out)->elts[1]), ((Params*)params)->x,
                    ((Params*)params)->y)))
    {
        free_vector_d(*grad_out); *grad_out = NULL;
        return err;
    }

    return NULL;
}

Error* langevin_accept_sample(const Vector_d* sample, void* params)
{
    double V;
    Error* err;

    set_params_from_pt(sample, params);
    if ((err = energy_func(&V, params)))
    {
        return err;
    }

    fprintf(stdout, "%7.4f  %7.4f  %8e\n",
            ((Params*)params)->x, ((Params*)params)->y, V);

    return NULL;
}

Error* hyper_accept_sample(const Vector_d* sample, void* params, double Vb)
{
    double V;
    Error* err;

    set_params_from_pt(sample, params);
    if ((err = energy_func(&V, params)))
    {
        return err;
    }

    fprintf(stdout, "%7.4f  %7.4f  %8e  %8e  %8e\n",
            ((Params*)params)->x, ((Params*)params)->y, V, Vb, Vb+V);

    return NULL;
}

Error* stochastic_dynamics_accept_sample
(
    const Vector_d* sample, 
    const Vector_d* momenta, 
    void*           params
)
{
    double E;
    double K;
    Error* err;

    set_params_from_pt(sample, params);
    if ((err = energy_func(&E, params)))
    {
        return err;
    }

    calc_vector_dot_prod_d(&K, momenta, momenta);
    K *= 0.5;

    fprintf(stdout, "%7.4f  %7.4f  %8e  %8e  %8e\n",
            ((Params*)params)->x, ((Params*)params)->y, E, K, E+K);

    return NULL;
}

Error* hybrid_mc_accept_sample
(
    const Vector_d* sample, 
    const Vector_d* momenta, 
    void*           params
)
{
    double E;
    double K;
    Error* err;

    set_params_from_pt(sample, params);
    if ((err = energy_func(&E, params)))
    {
        return err;
    }

    calc_vector_dot_prod_d(&K, momenta, momenta);
    K *= 0.5;

    fprintf(stdout, "%7.4f  %7.4f  %8e  %8e  %8e\n",
            ((Params*)params)->x, ((Params*)params)->y, E, K, E+K);

    return NULL;
}