|
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 00046 #include <config.h> 00047 00048 #include <stdlib.h> 00049 #include <stdio.h> 00050 #include <inttypes.h> 00051 #include <math.h> 00052 #include <assert.h> 00053 00054 #include <jwsc/base/error.h> 00055 #include <jwsc/matrix/matrix.h> 00056 00057 #include "potential.h" 00058 00059 00060 #define ERR_BUF_LEN 256 00061 00062 00064 static const float A[4] = {-200.0f, -100.0f, -170.0f, 15.0f}; 00065 00066 00068 static const float a[4] = {-1.0f, -1.0f, -6.5f, 0.7f}; 00069 00070 00072 static const float b[4] = {0.0f, 0.0f, 11.0f, 0.6f}; 00073 00074 00076 static const float c[4] = {-10.0f, -10.0f, -6.5f, 0.7f}; 00077 00078 00080 static const float X[4] = {1.0f, 0.0f, -0.5f, -1.0f}; 00081 00082 00084 static const float Y[4] = {0.0f, 0.5f, 1.5f, 1.0f}; 00085 00086 00088 static char err_buf[ERR_BUF_LEN] = {0}; 00089 00090 00106 Error* mullers_potential_f(float* V_out, float x, float y) 00107 { 00108 static char err_buf[256] = {0}; 00109 00110 uint8_t i; 00111 float V; 00112 00113 V = 0; 00114 00115 for (i = 0; i < 4; i++) 00116 { 00117 V += A[i] * exp(a[i] * (x - X[i]) * (x - X[i]) + 00118 b[i] * (x - X[i]) * (y - Y[i]) + 00119 c[i] * (y - Y[i]) * (y - Y[i])); 00120 } 00121 00122 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN 00123 if (isinf(V) || isnan(V)) 00124 { 00125 snprintf(err_buf, ERR_BUF_LEN-1, 00126 "Muller's potential precision error: x=%e y=%e", x, y); 00127 return JWSC_ECALC(err_buf); 00128 } 00129 #endif 00130 00131 *V_out = V; 00132 00133 return NULL; 00134 } 00135 00144 Error* mullers_potential_d(double* V_out, double x, double y) 00145 { 00146 uint8_t i; 00147 double V; 00148 00149 V = 0; 00150 00151 for (i = 0; i < 4; i++) 00152 { 00153 V += (double)A[i] * 00154 exp((double)a[i] * (x - (double)X[i]) * (x - (double)X[i]) + 00155 (double)b[i] * (x - (double)X[i]) * (y - (double)Y[i]) + 00156 (double)c[i] * (y - (double)Y[i]) * (y - (double)Y[i])); 00157 } 00158 00159 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN 00160 if (isinf(V) || isnan(V)) 00161 { 00162 snprintf(err_buf, ERR_BUF_LEN-1, 00163 "Muller's potential precision error: x=%e y=%e", x, y); 00164 return JWSC_ECALC(err_buf); 00165 } 00166 #endif 00167 00168 *V_out = V; 00169 00170 return NULL; 00171 } 00172 00198 Error* grad_mullers_potential_f 00199 ( 00200 float* dx_V_out, 00201 float* dy_V_out, 00202 float x, 00203 float y 00204 ) 00205 { 00206 uint8_t i; 00207 float dx_V, dy_V; 00208 00209 dx_V = 0; 00210 00211 for (i = 0; i < 4; i++) 00212 { 00213 dx_V += A[i] * (2.0 * a[i] * (x - X[i]) + b[i] * (y - Y[i])) * 00214 exp(a[i] * (x - X[i]) * (x - X[i]) + 00215 b[i] * (x - X[i]) * (y - Y[i]) + 00216 c[i] * (y - Y[i]) * (y - Y[i])); 00217 } 00218 00219 dy_V = 0; 00220 00221 for (i = 0; i < 4; i++) 00222 { 00223 dy_V += A[i] * (2.0 * c[i] * (y - Y[i]) + b[i] * (x - X[i])) * 00224 exp(a[i] * (x - X[i]) * (x - X[i]) + 00225 b[i] * (x - X[i]) * (y - Y[i]) + 00226 c[i] * (y - Y[i]) * (y - Y[i])); 00227 } 00228 00229 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN 00230 if (isinf(dx_V) || isnan(dx_V) || isinf(dy_V) || isnan(dy_V)) 00231 { 00232 snprintf(err_buf, ERR_BUF_LEN-1, 00233 "Gradient of Muller's potential precision error: x=%e y=%e", 00234 x, y); 00235 return JWSC_ECALC(err_buf); 00236 } 00237 #endif 00238 00239 *dx_V_out = dx_V; 00240 *dy_V_out = dy_V; 00241 00242 return NULL; 00243 } 00244 00258 Error* grad_mullers_potential_d 00259 ( 00260 double* dx_V_out, 00261 double* dy_V_out, 00262 double x, 00263 double y 00264 ) 00265 { 00266 uint8_t i; 00267 double dx_V, dy_V; 00268 00269 dx_V = 0; 00270 00271 for (i = 0; i < 4; i++) 00272 { 00273 dx_V += (double)A[i] * 00274 (2.0 * (double)a[i] * (x - (double)X[i]) + 00275 (double)b[i] * (y - (double)Y[i])) * 00276 exp((double)a[i] * (x - (double)X[i]) * (x - (double)X[i]) + 00277 (double)b[i] * (x - (double)X[i]) * (y - (double)Y[i]) + 00278 (double)c[i] * (y - (double)Y[i]) * (y - (double)Y[i])); 00279 } 00280 00281 dy_V = 0; 00282 00283 for (i = 0; i < 4; i++) 00284 { 00285 dy_V += (double)A[i] * 00286 (2.0 * (double)c[i] * (y - (double)Y[i]) + 00287 (double)b[i] * (x - (double)X[i])) * 00288 exp((double)a[i] * (x - (double)X[i]) * (x - (double)X[i]) + 00289 (double)b[i] * (x - (double)X[i]) * (y - (double)Y[i]) + 00290 (double)c[i] * (y - (double)Y[i]) * (y - (double)Y[i])); 00291 } 00292 00293 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN 00294 if (isinf(dx_V) || isnan(dx_V) || isinf(dy_V) || isnan(dy_V)) 00295 { 00296 snprintf(err_buf, ERR_BUF_LEN-1, 00297 "Gradient of Muller's potential precision error: x=%e y=%e", 00298 x, y); 00299 return JWSC_ECALC(err_buf); 00300 } 00301 #endif 00302 00303 *dx_V_out = dx_V; 00304 *dy_V_out = dy_V; 00305 00306 return NULL; 00307 } 00308 00332 Error* grad_matrix_mullers_potential_d 00333 ( 00334 Matrix_d** gg_out, 00335 double x, 00336 double y 00337 ) 00338 { 00339 double dx_V, dy_V; 00340 00341 Error* err; 00342 00343 create_matrix_d(gg_out, 2, 2); 00344 00345 if ((err = grad_mullers_potential_d(&dx_V, &dy_V, x, y))) 00346 { 00347 return err; 00348 } 00349 00350 (*gg_out)->elts[0][0] = dx_V*dx_V; 00351 (*gg_out)->elts[1][1] = dy_V*dy_V; 00352 (*gg_out)->elts[0][1] = (*gg_out)->elts[1][0] = dx_V*dy_V; 00353 00354 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN 00355 if (isinf(dx_V*dx_V) || isnan(dx_V*dx_V) || 00356 isinf(dy_V*dy_V) || isnan(dy_V*dy_V) || 00357 isinf(dx_V*dy_V) || isnan(dx_V*dy_V)) 00358 { 00359 snprintf(err_buf, ERR_BUF_LEN-1, 00360 "Gradient of Muller's potential matrix precision error: x=%e y=%e", 00361 x, y); 00362 return JWSC_ECALC(err_buf); 00363 } 00364 #endif 00365 00366 return NULL; 00367 } 00368 00392 Error* hessian_mullers_potential_f 00393 ( 00394 Matrix_f** H_V_out, 00395 float x, 00396 float y 00397 ) 00398 { 00399 uint8_t i; 00400 float* dx2_V; 00401 float* dy2_V; 00402 float* dxy_V; 00403 00404 Matrix_f* H_V; 00405 00406 create_zero_matrix_f(H_V_out, 2, 2); 00407 H_V = *H_V_out; 00408 00409 dx2_V = &(H_V->elts[0][0]); 00410 00411 for (i = 0; i < 4; i++) 00412 { 00413 *dx2_V += A[i] * 00414 (2.0*a[i] + (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i])) * 00415 (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) * 00416 exp(a[i] * (x - X[i]) * (x - X[i]) + 00417 b[i] * (x - X[i]) * (y - Y[i]) + 00418 c[i] * (y - Y[i]) * (y - Y[i])); 00419 } 00420 00421 dy2_V = &(H_V->elts[1][1]); 00422 00423 for (i = 0; i < 4; i++) 00424 { 00425 *dy2_V += A[i] * 00426 (2.0*c[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) * 00427 (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i]))) * 00428 exp(a[i] * (x - X[i]) * (x - X[i]) + 00429 b[i] * (x - X[i]) * (y - Y[i]) + 00430 c[i] * (y - Y[i]) * (y - Y[i])); 00431 } 00432 00433 dxy_V = &(H_V->elts[0][1]); 00434 00435 for (i = 0; i < 4; i++) 00436 { 00437 *dxy_V += A[i] * 00438 (b[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) * 00439 (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) * 00440 exp(a[i] * (x - X[i]) * (x - X[i]) + 00441 b[i] * (x - X[i]) * (y - Y[i]) + 00442 c[i] * (y - Y[i]) * (y - Y[i])); 00443 } 00444 00445 H_V->elts[1][0] = H_V->elts[0][1]; 00446 00447 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN 00448 if (isinf(*dx2_V) || isnan(*dx2_V) || isinf(*dy2_V) || isnan(*dy2_V) || 00449 isinf(*dxy_V) || isnan(*dxy_V)) 00450 { 00451 snprintf(err_buf, ERR_BUF_LEN-1, 00452 "Hessian of Muller's potential precision error: x=%e y=%e", 00453 x, y); 00454 return JWSC_ECALC(err_buf); 00455 } 00456 #endif 00457 00458 return NULL; 00459 } 00460 00472 Error* hessian_mullers_potential_d 00473 ( 00474 Matrix_d** H_V_out, 00475 double x, 00476 double y 00477 ) 00478 { 00479 uint8_t i; 00480 double* dx2_V; 00481 double* dy2_V; 00482 double* dxy_V; 00483 00484 Matrix_d* H_V; 00485 00486 create_zero_matrix_d(H_V_out, 2, 2); 00487 H_V = *H_V_out; 00488 00489 dx2_V = &(H_V->elts[0][0]); 00490 00491 for (i = 0; i < 4; i++) 00492 { 00493 *dx2_V += A[i] * 00494 (2.0*a[i] + (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i])) * 00495 (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) * 00496 exp(a[i] * (x - X[i]) * (x - X[i]) + 00497 b[i] * (x - X[i]) * (y - Y[i]) + 00498 c[i] * (y - Y[i]) * (y - Y[i])); 00499 } 00500 00501 dy2_V = &(H_V->elts[1][1]); 00502 00503 for (i = 0; i < 4; i++) 00504 { 00505 *dy2_V += A[i] * 00506 (2.0*c[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) * 00507 (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i]))) * 00508 exp(a[i] * (x - X[i]) * (x - X[i]) + 00509 b[i] * (x - X[i]) * (y - Y[i]) + 00510 c[i] * (y - Y[i]) * (y - Y[i])); 00511 } 00512 00513 dxy_V = &(H_V->elts[0][1]); 00514 00515 for (i = 0; i < 4; i++) 00516 { 00517 *dxy_V += A[i] * 00518 (b[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) * 00519 (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) * 00520 exp(a[i] * (x - X[i]) * (x - X[i]) + 00521 b[i] * (x - X[i]) * (y - Y[i]) + 00522 c[i] * (y - Y[i]) * (y - Y[i])); 00523 } 00524 00525 H_V->elts[1][0] = H_V->elts[0][1]; 00526 00527 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN 00528 if (isinf(*dx2_V) || isnan(*dx2_V) || isinf(*dy2_V) || isnan(*dy2_V) || 00529 isinf(*dxy_V) || isnan(*dxy_V)) 00530 { 00531 snprintf(err_buf, ERR_BUF_LEN-1, 00532 "Hessian of Muller's potential precision error: x=%e y=%e", 00533 x, y); 00534 return JWSC_ECALC(err_buf); 00535 } 00536 #endif 00537 00538 return NULL; 00539 } 00540
1.7.4