Muller
sample Muller's potential
sampler.c
Go to the documentation of this file.
00001 /*
00002  * This work is licensed under a Creative Commons 
00003  * Attribution-Noncommercial-Share Alike 3.0 United States License.
00004  * 
00005  *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
00006  * 
00007  * You are free:
00008  * 
00009  *    to Share - to copy, distribute, display, and perform the work
00010  *    to Remix - to make derivative works
00011  * 
00012  * Under the following conditions:
00013  * 
00014  *    Attribution. You must attribute the work in the manner specified by the
00015  *    author or licensor (but not in any way that suggests that they endorse you
00016  *    or your use of the work).
00017  * 
00018  *    Noncommercial. You may not use this work for commercial purposes.
00019  * 
00020  *    Share Alike. If you alter, transform, or build upon this work, you may
00021  *    distribute the resulting work only under the same or similar license to
00022  *    this one.
00023  * 
00024  * For any reuse or distribution, you must make clear to others the license
00025  * terms of this work. The best way to do this is by including this header.
00026  * 
00027  * Any of the above conditions can be waived if you get permission from the
00028  * copyright holder.
00029  * 
00030  * Apart from the remix rights granted under this license, nothing in this
00031  * license impairs or restricts the author's moral rights.
00032  */
00033 
00034 
00047 #include <config.h>
00048 
00049 #include <stdlib.h>
00050 #include <stdio.h>
00051 #include <math.h>
00052 #include <assert.h>
00053 #include <inttypes.h>
00054 
00055 #include <jwsc/base/error.h>
00056 #include <jwsc/base/limits.h>
00057 #include <jwsc/prob/pdf.h>
00058 #include <jwsc/vector/vector.h>
00059 #include <jwsc/vector/vector_math.h>
00060 #include <jwsc/matrix/matrix.h>
00061 #include <jwsc/matrix/matrix_math.h>
00062 
00063 #include "potential.h"
00064 #include "sampler.h"
00065 
00066 
00079 Error* metropolis_hastings
00080 (
00081     uint32_t iterations,
00082     double   start_x,
00083     double   start_y,
00084     double   x_sigma,
00085     double   x_min,
00086     double   x_max,
00087     double   y_sigma,
00088     double   y_min,
00089     double   y_max
00090 )
00091 {
00092     uint32_t i;
00093     uint8_t  accepted;
00094 
00095     double log_alpha;
00096     double log_u;
00097 
00098     double x_accepted;
00099     double y_accepted;
00100     double V_accepted;
00101 
00102     double x_proposed;
00103     double y_proposed;
00104     double V_proposed;
00105 
00106     double q_accepted_lp;
00107     double q_proposed_lp;
00108 
00109     Error* err;
00110 
00111     x_accepted = start_x;
00112     y_accepted = start_y;
00113     x_proposed = x_accepted;
00114     y_proposed = y_accepted;
00115 
00116     if ((err = mullers_potential_d(&V_accepted, x_accepted, y_accepted)))
00117     {
00118         return err;
00119     }
00120     V_proposed = V_accepted;
00121 
00122     accepted = 1;
00123 
00124     for (i = 0; i < iterations; i++)
00125     {
00126         fprintf(stdout, "%7.4f  %7.4f  %8e  %7.4f  %7.4f  %8e  %1u\n",
00127                 x_accepted,
00128                 y_accepted,
00129                 V_accepted, 
00130                 x_proposed,
00131                 y_proposed,
00132                 V_proposed,
00133                 accepted);
00134 
00135         x_proposed = sample_gaussian_pdf_d(x_accepted, x_sigma, x_min, x_max);
00136         y_proposed = sample_gaussian_pdf_d(y_accepted, y_sigma, y_min, y_max);
00137 
00138         if ((err = mullers_potential_d(&V_proposed, x_proposed, y_proposed)))
00139         {
00140             return err;
00141         }
00142 
00143         q_proposed_lp  = 0;
00144         q_proposed_lp += gaussian_pdf_d(x_accepted, x_sigma, x_proposed);
00145         q_proposed_lp += gaussian_pdf_d(y_accepted, y_sigma, y_proposed);
00146 
00147         q_accepted_lp  = 0;
00148         q_accepted_lp += gaussian_pdf_d(x_proposed, x_sigma, x_accepted);
00149         q_accepted_lp += gaussian_pdf_d(y_proposed, y_sigma, y_accepted);
00150 
00151         log_alpha  = 0;
00152         log_alpha += log((-V_proposed < JWSC_MIN_LOG_ARG) ? 
00153                          JWSC_MIN_LOG_ARG : -V_proposed);
00154         log_alpha -= log((-V_accepted < JWSC_MIN_LOG_ARG) ?
00155                          JWSC_MIN_LOG_ARG : -V_accepted);
00156         log_alpha += q_accepted_lp;
00157         log_alpha -= q_proposed_lp;
00158 
00159         log_u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1.0));
00160 
00161         accepted = 0;
00162 
00163         if (log_u <= log_alpha)
00164         {
00165             accepted = 1;
00166             x_accepted = x_proposed;
00167             y_accepted = y_proposed;
00168             V_accepted = V_proposed;
00169         }
00170     }
00171 
00172     return NULL;
00173 }
00174 
00175 
00176 void get_pt_from_params(Vector_d** p_out, const void* params)
00177 {
00178     create_vector_d(p_out, 2);
00179     (*p_out)->elts[0] = ((Params*)params)->x;
00180     (*p_out)->elts[1] = ((Params*)params)->y;
00181 }
00182 
00183 Error* set_params_from_pt(const Vector_d* p, void* params)
00184 {
00185     ((Params*)params)->x = p->elts[0];
00186     ((Params*)params)->y = p->elts[1];
00187 
00188     return NULL;
00189 }
00190 
00191 Error* energy_func(double* energy_out, void* params)
00192 {
00193     Error* err;
00194 
00195     if ((err = mullers_potential_d(energy_out, ((Params*)params)->x,
00196                     ((Params*)params)->y)))
00197     {
00198         return err;
00199     }
00200 
00201     return NULL;
00202 }
00203 
00204 Error* grad_energy_func(Vector_d** grad_out, void* params)
00205 {
00206     Error* err;
00207 
00208     create_vector_d(grad_out, 2);
00209 
00210     if ((err = grad_mullers_potential_d(&((*grad_out)->elts[0]), 
00211                     &((*grad_out)->elts[1]), ((Params*)params)->x,
00212                     ((Params*)params)->y)))
00213     {
00214         free_vector_d(*grad_out); *grad_out = NULL;
00215         return err;
00216     }
00217 
00218     return NULL;
00219 }
00220 
00221 Error* langevin_accept_sample(const Vector_d* sample, void* params)
00222 {
00223     double V;
00224     Error* err;
00225 
00226     set_params_from_pt(sample, params);
00227     if ((err = energy_func(&V, params)))
00228     {
00229         return err;
00230     }
00231 
00232     fprintf(stdout, "%7.4f  %7.4f  %8e\n",
00233             ((Params*)params)->x, ((Params*)params)->y, V);
00234 
00235     return NULL;
00236 }
00237 
00238 Error* hyper_accept_sample(const Vector_d* sample, void* params, double Vb)
00239 {
00240     double V;
00241     Error* err;
00242 
00243     set_params_from_pt(sample, params);
00244     if ((err = energy_func(&V, params)))
00245     {
00246         return err;
00247     }
00248 
00249     fprintf(stdout, "%7.4f  %7.4f  %8e  %8e  %8e\n",
00250             ((Params*)params)->x, ((Params*)params)->y, V, Vb, Vb+V);
00251 
00252     return NULL;
00253 }
00254 
00255 Error* stochastic_dynamics_accept_sample
00256 (
00257     const Vector_d* sample, 
00258     const Vector_d* momenta, 
00259     void*           params
00260 )
00261 {
00262     double E;
00263     double K;
00264     Error* err;
00265 
00266     set_params_from_pt(sample, params);
00267     if ((err = energy_func(&E, params)))
00268     {
00269         return err;
00270     }
00271 
00272     calc_vector_dot_prod_d(&K, momenta, momenta);
00273     K *= 0.5;
00274 
00275     fprintf(stdout, "%7.4f  %7.4f  %8e  %8e  %8e\n",
00276             ((Params*)params)->x, ((Params*)params)->y, E, K, E+K);
00277 
00278     return NULL;
00279 }
00280 
00281 Error* hybrid_mc_accept_sample
00282 (
00283     const Vector_d* sample, 
00284     const Vector_d* momenta, 
00285     void*           params
00286 )
00287 {
00288     double E;
00289     double K;
00290     Error* err;
00291 
00292     set_params_from_pt(sample, params);
00293     if ((err = energy_func(&E, params)))
00294     {
00295         return err;
00296     }
00297 
00298     calc_vector_dot_prod_d(&K, momenta, momenta);
00299     K *= 0.5;
00300 
00301     fprintf(stdout, "%7.4f  %7.4f  %8e  %8e  %8e\n",
00302             ((Params*)params)->x, ((Params*)params)->y, E, K, E+K);
00303 
00304     return NULL;
00305 }