diff --git a/msmtools/estimation/dense/sample_rev.c b/msmtools/estimation/dense/sample_rev.c index f0901b36..7d89be80 100644 --- a/msmtools/estimation/dense/sample_rev.c +++ b/msmtools/estimation/dense/sample_rev.c @@ -32,38 +32,16 @@ #include "_ranlib.h" #include "sample_rev.h" +#include "util.h" #define _square(x) x*x -/** - Helper function, tests if x is numerically positive - - :param x: - :return: -*/ -int _is_positive(double x) -{ - double eps = 1e-8; -#if _MSC_VER && !__INTEL_COMPILER - if(x >=eps && ! (!_finite(x) && !_isnan(x))) -#else - if (x >= eps && !isinf(x) && !isnan(x)) -#endif - return 1; - else - return 0; -} - int _accept_step(double log_prob_old, double log_prob_new) { if (log_prob_new > log_prob_old) // this is faster return 1; -#ifdef _MSC_VER && !__INTEL_COMPILER - if (genunf(0,1) < exp( __min( log_prob_new-log_prob_old, 0 ) )) -#else - if (genunf(0,1) < exp( fmin( log_prob_new-log_prob_old, 0 ) )) -#endif + if (genunf(0,1) < exp( my_fmin( log_prob_new-log_prob_old, 0 ) )) return 1; return 0; } @@ -103,11 +81,11 @@ double _update_step(double v0, double v1, double v2, double c0, double c1, doubl // about 1.5 sec: gamma and normf generation // about 1 sec: logs+exps in else blocks - if (_is_positive(k) && _is_positive(theta)) + if (is_positive(k) && is_positive(theta)) { v0_new = gengam(1.0/theta,k); log_v0_new = log(v0_new); - if (_is_positive(v0_new)) + if (is_positive(v0_new)) { if (v0 == 0) { @@ -121,7 +99,7 @@ double _update_step(double v0, double v1, double v2, double c0, double c1, doubl log_prob_old = (c0-1) * log_v0 - c1 * log(v0+v1) - c2 * log(v0+v2); log_prob_old -= (k-1) * log_v0 - v0/theta; if (_accept_step(log_prob_old, log_prob_new)) - //if (genunf(0,1) < exp( fmin( log_prob_new-log_prob_old, 0 ) )) + //if (genunf(0,1) < exp( my_fmin( log_prob_new-log_prob_old, 0 ) )) { v0 = v0_new; log_v0 = log_v0_new; @@ -132,7 +110,7 @@ double _update_step(double v0, double v1, double v2, double c0, double c1, doubl v0_new = v0 * exp( random_walk_stepsize * snorm() ); log_v0_new = log(v0_new); - if (_is_positive(v0_new)) + if (is_positive(v0_new)) { if (v0 == 0) { @@ -144,7 +122,7 @@ double _update_step(double v0, double v1, double v2, double c0, double c1, doubl log_prob_new = c0 * log_v0_new - c1 * log(v0_new + v1) - c2 * log(v0_new + v2); log_prob_old = c0 * log_v0 - c1 * log(v0 + v1) - c2 * log(v0 + v2); if (_accept_step(log_prob_old, log_prob_new)) - //if (genunf(0,1) < exp( fmin( log_prob_new-log_prob_old, 0 ) )) + //if (genunf(0,1) < exp( my_fmin( log_prob_new-log_prob_old, 0 ) )) { v0 = v0_new; log_v0 = log_v0_new; @@ -262,11 +240,11 @@ void _update(double* C, double* sumC, double* X, int n, int n_step) { if (i == j) { - if (_is_positive(C[i*n+i]) && _is_positive(sumC[i] - C[i*n+i])) + if (is_positive(C[i*n+i]) && is_positive(sumC[i] - C[i*n+i])) { tmp1 = genbet(C[i*n+i], sumC[i] - C[i*n+i]); tmp2 = tmp1 / (1-tmp1) * (_sum_row(X, n, i) - X[i*n+i]); - if (_is_positive(tmp2)) + if (is_positive(tmp2)) { X[i*n+i] = tmp2; } @@ -363,11 +341,11 @@ void _update_sparse(double* C, double* sumC, double* X, double* sumX, int* I, in j = J[k]; if (i == j) { - if (_is_positive(C[i*n+i]) && _is_positive(sumC[i] - C[i*n+i])) + if (is_positive(C[i*n+i]) && is_positive(sumC[i] - C[i*n+i])) { tmp1 = genbet(C[i*n+i], sumC[i] - C[i*n+i]); tmp2 = tmp1 / (1-tmp1) * (sumX[i] - X[i*n+i]); - if (_is_positive(tmp2)) + if (is_positive(tmp2)) { sumX[i] += tmp2 - X[i*n+i]; // update sumX X[i*n+i] = tmp2; diff --git a/msmtools/estimation/dense/sample_rev.h b/msmtools/estimation/dense/sample_rev.h index 02c80a68..e8d743e7 100644 --- a/msmtools/estimation/dense/sample_rev.h +++ b/msmtools/estimation/dense/sample_rev.h @@ -26,7 +26,9 @@ /* * moduleauthor:: F. Noe */ -int _is_positive(double x); +#ifndef _sample_rev_h_ +#define _sample_rev_h_ + double _update_step(double v0, double v1, double v2, double c0, double c1, double c2, int random_walk_stepsize); double _sum_all(double* X, int n); void _normalize_all(double* X, int n); @@ -43,3 +45,6 @@ void _update_sparse_speedtest(double* C, double* sumC, double* X, double* sumX, void _update_sparse(double* C, double* sumC, double* X, double* sumX, int* I, int* J, int n, int n_idx, int n_step); void _generate_row_indexes(int* I, int n, int n_idx, int* row_indexes); void _normalize_all(double* X, int n); + + +#endif diff --git a/msmtools/estimation/dense/sample_revpi.c b/msmtools/estimation/dense/sample_revpi.c index 58f9d9ac..98ae255c 100644 --- a/msmtools/estimation/dense/sample_revpi.c +++ b/msmtools/estimation/dense/sample_revpi.c @@ -33,19 +33,8 @@ #include "_ranlib.h" #include "sample_revpi.h" +#include "util.h" -int is_positive(double x) -{ - double eps = 1e-15; -#if _MSC_VER && !__INTEL_COMPILER - if(x >=eps && ! (!_finite(x) && !_isnan(x))) -#else - if (x >= eps && !isinf(x) && !isnan(x)) -#endif - return 1; - else - return 0; -} double f(double v, double s, double a1, double a2, double a3) @@ -170,11 +159,7 @@ sample_quad(double xkl, double xkk, double xll, // Metropolis step U = genunf(0.0, 1.0); -#if _MSC_VER && !__INTEL_COMPILER - if(log(U) < __min(0.0, q)) -#else - if(log(U) < fmin(0.0, q)) -#endif + if(log(U) < my_fmin(0.0, q)) { return s2*w/(1.0+w); } @@ -238,11 +223,7 @@ sample_quad_rw(double xkl, double xkk, double xll, q = qacc_rw(w, v, s, a1, a2, a3); //Metropolis step U = genunf(0.0, 1.0); -#if _MSC_VER && !__INTEL_COMPILER - if (log(U) < __min(0.0, q)) -#else - if(log(U) < fmin(0.0, q)) -#endif + if(log(U) < my_fmin(0.0, q)) { return s2*w/(1.0+w); } diff --git a/msmtools/estimation/dense/util.h b/msmtools/estimation/dense/util.h new file mode 100644 index 00000000..9b0baf7f --- /dev/null +++ b/msmtools/estimation/dense/util.h @@ -0,0 +1,50 @@ +#ifndef __util_h_ +#define __util_h_ + +#include +#include +#include +#include + + +int my_isinf(double x) { +#if _MSC_VER && !__INTEL_COMPILER + return ! _finite(x); +#else + return isinf(x); +#endif +} + +int my_isnan(double x) { +#if _MSC_VER && !__INTEL_COMPILER + return _isnan(x); +#else + return isnan(x); +#endif +} + +/** + Helper function, tests if x is numerically positive + + :param x: + :return: +*/ +int is_positive(double x) +{ + double eps = 1e-8; + if (x >= eps && !my_isinf(x) && !my_isnan(x)) + return 1; + else + return 0; +} + +double my_fmin(double a, double b) { +#if _MSC_VER && !__INTEL_COMPILER + return __min(a,b); +#else + return fmin(a,b); +#endif +} + + +#endif