Skip to content

Commit

Permalink
Merge pull request #42 from rcppsmc/feature/more_std
Browse files Browse the repository at this point in the history
use std::pow instead of pow to fix solaris build
  • Loading branch information
eddelbuettel committed Mar 18, 2018
2 parents 14b391e + 659ad6b commit 3211f50
Show file tree
Hide file tree
Showing 9 changed files with 28 additions and 20 deletions.
8 changes: 8 additions & 0 deletions ChangeLog
@@ -1,5 +1,13 @@
2018-03-17 Dirk Eddelbuettel <edd@debian.org>

* inst/include/history.h: Use std::pow
* inst/include/sampler.h: Idem
* src/LinReg.cpp: Idem
* src/LinReg_LA_adapt.cpp: Idem
* src/nonLinPMMH.cpp: Idem
* src/pflineart.cpp: Idem
* src/pfnonlinbs.cpp: Idem

* DESCRIPTION (Version, Date): New minor version

2018-03-16 Adam M. Johansen <adam.johansen@gmail.com>
Expand Down
2 changes: 1 addition & 1 deletion inst/include/history.h
Expand Up @@ -168,7 +168,7 @@ namespace smc {
long double wSum = expl(stableLogSumWeights(pop.GetLogWeight()));
for(long i =0; i < number; i++)
{
rValue += expl(pop.GetLogWeightN(i)) * pow(static_cast<long double>(pIntegrand(lTime, pop.GetValueN(i), pAuxiliary) - Expectation),2.0);
rValue += expl(pop.GetLogWeightN(i)) * std::pow(static_cast<long double>(pIntegrand(lTime, pop.GetValueN(i), pAuxiliary) - Expectation),2.0);
}

rValue /= wSum;
Expand Down
2 changes: 1 addition & 1 deletion inst/include/sampler.h
Expand Up @@ -508,7 +508,7 @@ namespace smc {
current_expt = it->Integrate(lTime, pIntegrand, pAuxiliary);
current_var = it->Integrate_Var(lTime, pIntegrand, current_expt, pAuxiliary);
width = static_cast<long double>(pWidth(lTime,pAuxiliary));
rValue += width/2.0 * (previous_expt + current_expt) - pow(width,2.0)/12.0*(current_var - previous_var);
rValue += width/2.0 * (previous_expt + current_expt) - std::pow(width,2.0)/12.0*(current_var - previous_var);
lTime++;
previous_expt = current_expt;
previous_var = current_var;
Expand Down
14 changes: 7 additions & 7 deletions src/LinReg.cpp
Expand Up @@ -35,7 +35,7 @@ namespace LinReg {
arma::mat covRW("2500 -2.5 0.03; -2.5 130.0 0.0; 0.03 0.0 0.04");
arma::mat cholCovRW = arma::chol(covRW);
const double a_prior = 3.0;
const double b_prior = pow(2.0*300.0*300.0,-1.0);
const double b_prior = std::pow(2.0*300.0*300.0,-1.0);
}

using namespace std;
Expand Down Expand Up @@ -92,8 +92,8 @@ namespace LinReg {
double logWeight(long lTime, const rad_state & value){

double mean_reg = value.theta(0) + value.theta(1)*(data.x(lTime) - mean_x);
double sigma = pow(expl(value.theta(2)),0.5);
return -log(sigma) - pow(data.y(lTime) - mean_reg,2.0)/(2.0*sigma*sigma) -0.5*log(2.0*M_PI);
double sigma = std::pow(expl(value.theta(2)),0.5);
return -log(sigma) - std::pow(data.y(lTime) - mean_reg,2.0)/(2.0*sigma*sigma) -0.5*log(2.0*M_PI);

}

Expand All @@ -103,15 +103,15 @@ namespace LinReg {
/// \param value The state to consider
double logPosterior(long lTime, const rad_state & value){

double log_prior = -log(1000.0)- pow(value.theta(0) - 3000.0,2.0)/(2.0*1000.0*1000.0) -log(100.0)- pow(value.theta(1) - 185.0,2.0)/(2.0*100.0*100.0) + value.theta(2)-1.0/b_prior/expl(value.theta(2)) -value.theta(2)*(a_prior+1.0);
double log_prior = -log(1000.0)- std::pow(value.theta(0) - 3000.0,2.0)/(2.0*1000.0*1000.0) -log(100.0)- std::pow(value.theta(1) - 185.0,2.0)/(2.0*100.0*100.0) + value.theta(2)-1.0/b_prior/expl(value.theta(2)) -value.theta(2)*(a_prior+1.0);

double sigma = pow(expl(value.theta(2)),0.5);
double sigma = std::pow(expl(value.theta(2)),0.5);

double log_normpdf;

if (lTime==0){
double mean_reg = value.theta(0) + value.theta(1)*(data.x(0) - mean_x);
log_normpdf = -log(sigma) - pow(data.y(0) - mean_reg,2.0)/(2.0*sigma*sigma) -0.5*log(2.0*M_PI);
log_normpdf = -log(sigma) - std::pow(data.y(0) - mean_reg,2.0)/(2.0*sigma*sigma) -0.5*log(2.0*M_PI);
} else{
arma::vec mean_reg = value.theta(0) + value.theta(1)*(data.x.rows(0,lTime) - mean_x);
log_normpdf = arma::sum(-log(sigma) - pow(data.y.rows(0,lTime) - mean_reg,2.0)/(2.0*sigma*sigma) -0.5*log(2.0*M_PI));
Expand All @@ -131,7 +131,7 @@ namespace LinReg {
// drawing from the prior
value.theta(0) = R::rnorm(3000.0,1000.0);
value.theta(1) = R::rnorm(185.0,100.0);
value.theta(2) = log(pow(R::rgamma(3,pow(2.0*300.0*300.0,-1.0)),-1.0));
value.theta(2) = log(std::pow(R::rgamma(3,std::pow(2.0*300.0*300.0,-1.0)),-1.0));

logweight = logWeight(0, value);
}
Expand Down
4 changes: 2 additions & 2 deletions src/LinReg_LA.cpp
Expand Up @@ -125,7 +125,7 @@ namespace LinReg_LA {
/// \param value The state to consider
double logLikelihood(const rad_state & value){

double sigma = pow(expl(value.theta(2)),0.5);
double sigma = std::pow(expl(value.theta(2)),0.5);
arma::vec mean_reg = value.theta(0) + value.theta(1)*(data.x - mean_x);
return arma::sum(-log(sigma) - pow(data.y - mean_reg,2.0)/(2.0*sigma*sigma) -0.5*log(2.0*M_PI));

Expand Down Expand Up @@ -185,4 +185,4 @@ namespace LinReg_LA {
}
return FALSE;
}
}
}
8 changes: 4 additions & 4 deletions src/LinReg_LA_adapt.cpp
Expand Up @@ -26,7 +26,7 @@

namespace LinReg_LA_adapt {
const double a_prior = 3.0;
const double b_prior = pow(2.0*300.0*300.0,-1.0);
const double b_prior = std::pow(2.0*300.0*300.0,-1.0);
}

using namespace std;
Expand Down Expand Up @@ -118,15 +118,15 @@ namespace LinReg_LA_adapt {
/// \param value The state to consider
double logLikelihood(const rad_state & value){

double sigma = pow(expl(value.theta(2)),0.5);
double sigma = std::pow(expl(value.theta(2)),0.5);
arma::vec mean_reg = value.theta(0) + value.theta(1)*(data.x - mean_x);
return arma::sum(-log(sigma) - pow(data.y - mean_reg,2.0)/(2.0*sigma*sigma) -0.5*log(2.0*M_PI));

}
///The function corresponding to the (unnormalised) log prior at a specified position
/// \param value The state to consider
double logPrior(const rad_state & value){
return -log(1000.0)- pow(value.theta(0) - 3000.0,2.0)/(2.0*1000.0*1000.0) -log(100.0)- pow(value.theta(1) - 185.0,2.0)/(2.0*100.0*100.0) + value.theta(2)-1.0/b_prior/expl(value.theta(2)) -value.theta(2)*(a_prior+1.0);
return -log(1000.0)- std::pow(value.theta(0) - 3000.0,2.0)/(2.0*1000.0*1000.0) -log(100.0)- std::pow(value.theta(1) - 185.0,2.0)/(2.0*100.0*100.0) + value.theta(2)-1.0/b_prior/expl(value.theta(2)) -value.theta(2)*(a_prior+1.0);
}

///A function to initialise a particle
Expand All @@ -140,7 +140,7 @@ namespace LinReg_LA_adapt {
value.theta = arma::zeros(3);
value.theta(0) = R::rnorm(3000.0,1000.0);
value.theta(1) = R::rnorm(185.0,100.0);
value.theta(2) = log(pow(R::rgamma(3,pow(2.0*300.0*300.0,-1.0)),-1.0));
value.theta(2) = log(std::pow(R::rgamma(3,pow(2.0*300.0*300.0,-1.0)),-1.0));
value.loglike = logLikelihood(value);
value.logprior = logPrior(value);
logweight = params.GetTemp(0)*value.loglike;
Expand Down
6 changes: 3 additions & 3 deletions src/nonLinPMMH.cpp
Expand Up @@ -129,7 +129,7 @@ namespace nonLinPMMH {
void fInitialise(double & X, double & logweight, smc::nullParams & param)
{
X = R::rnorm(0.0,sqrt(5.0));
double mean = pow(X,2)/20.0;
double mean = std::pow(X,2)/20.0;
logweight = R::dnorm(y(0),mean,theta_prop.sigw,TRUE);
}

Expand All @@ -141,8 +141,8 @@ namespace nonLinPMMH {
/// \param param Additional algorithm parameters
void fMove(long lTime, double & X, double & logweight, smc::nullParams & param)
{
X = X/2.0 + 25.0*X/(1+pow(X,2)) + 8*cos(1.2*(lTime+1)) + R::rnorm(0.0,theta_prop.sigv);
double mean = pow(X,2)/20.0;
X = X/2.0 + 25.0*X/(1+std::pow(X,2)) + 8*cos(1.2*(lTime+1)) + R::rnorm(0.0,theta_prop.sigv);
double mean = std::pow(X,2)/20.0;
logweight += R::dnorm(y(lTime),mean,theta_prop.sigw,TRUE);
}

Expand Down
2 changes: 1 addition & 1 deletion src/pflineart.cpp
Expand Up @@ -134,7 +134,7 @@ namespace pflineart {
/// \param X The state to consider
double logLikelihood(long lTime, const cv_state & X)
{
return - 0.5 * (nu_y + 1.0) * (log(1 + pow((X.x_pos - y.x_pos(lTime))/scale_y,2) / nu_y) + log(1 + pow((X.y_pos - y.y_pos(lTime))/scale_y,2) / nu_y));
return - 0.5 * (nu_y + 1.0) * (log(1 + std::pow((X.x_pos - y.x_pos(lTime))/scale_y,2) / nu_y) + log(1 + std::pow((X.y_pos - y.y_pos(lTime))/scale_y,2) / nu_y));
}

///A function to initialise particles
Expand Down
2 changes: 1 addition & 1 deletion src/pfnonlinbs.cpp
Expand Up @@ -83,7 +83,7 @@ namespace nonlinbs {
/// \param lTime The current time (i.e. the index of the current distribution)
/// \param X The state to consider
double logLikelihood(long lTime, const double & x) {
return -0.5 * pow(y(lTime) - x*x*scale_y,2.0) / var_y;
return -0.5 * std::pow(y(lTime) - x*x*scale_y,2.0) / var_y;
}

///A function to initialise particles
Expand Down

0 comments on commit 3211f50

Please sign in to comment.