Muller
sample Muller's potential
|
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 }