Skip to content

Commit

Permalink
Merge pull request #15 from sor16/Rcpp-to-speed-up-mcmc
Browse files Browse the repository at this point in the history
Rcpp to speed up mcmc (and various other changes)
  • Loading branch information
RafaelVias authored Sep 4, 2024
2 parents 064d754 + 61f44d9 commit dd7eb96
Show file tree
Hide file tree
Showing 342 changed files with 36,686 additions and 16,397 deletions.
Binary file modified .DS_Store
Binary file not shown.
14 changes: 9 additions & 5 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
Package: bdrc
Title: Bayesian Discharge Rating Curves
Version: 1.1.0.9000
Version: 2.0.0
Authors@R: c(person("Birgir", "Hrafnkelsson", email = "[email protected]", role = c("aut","cph"),comment=c(ORCID="0000-0003-1864-9652")),
person("Solvi", "Rognvaldsson", email = "[email protected]", role = c("aut","cre"),comment=c(ORCID="0000-0002-4376-3361")),
person("Axel Orn", "Jansson", email = "axelorn94@gmail.com", role = c("aut")),
person("Rafael Daníel", "Vias", email = "raffidv@gmail.com", role = c("aut"),comment=c(ORCID="0009-0007-2601-6800"))
person("Solvi", "Rognvaldsson", email = "[email protected]", role = c("aut"),comment=c(ORCID="0000-0002-4376-3361")),
person(given = "Rafael Daníel", family = "Vias", email = "raffidv@gmail.com", role = c("aut", "cre"), comment=c(ORCID = "0009-0007-2601-6800")),
person("Axel Orn", "Jansson", email = "axelorn94@gmail.com", role = c("aut"))
)
Maintainer: Solvi Rognvaldsson <solviro@gmail.com>
Maintainer: Rafael Daníel Vias <raffidv@gmail.com>
Description: Fits a discharge rating curve based on the power-law and the generalized power-law from data on paired stage and discharge measurements in a given river using a Bayesian hierarchical model as described in Hrafnkelsson et al. (2020) <arXiv:2010.04769>.
Depends: R (>= 3.5.0)
License: MIT + file LICENSE
Expand All @@ -16,6 +16,7 @@ Imports:
ggplot2,
grid,
gridExtra,
Rcpp,
rlang,
scales
Suggests:
Expand All @@ -28,3 +29,6 @@ VignetteBuilder: knitr
URL: https://sor16.github.io/bdrc
BugReports: https://github.com/sor16/bdrc/issues
Encoding: UTF-8
LinkingTo:
Rcpp,
RcppArmadillo
4 changes: 4 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ export(plm)
export(plm0)
export(spread_draws)
export(tournament)
importFrom(Rcpp,sourceCpp)
importFrom(ggplot2,"%+replace%")
importFrom(ggplot2,aes)
importFrom(ggplot2,autoplot)
Expand Down Expand Up @@ -71,6 +72,7 @@ importFrom(ggplot2,ggtitle)
importFrom(ggplot2,guide_legend)
importFrom(ggplot2,guides)
importFrom(ggplot2,label_parsed)
importFrom(ggplot2,labs)
importFrom(ggplot2,margin)
importFrom(ggplot2,scale_color_manual)
importFrom(ggplot2,scale_colour_manual)
Expand Down Expand Up @@ -111,3 +113,5 @@ importFrom(stats,rnorm)
importFrom(stats,runif)
importFrom(stats,var)
importFrom(utils,askYesNo)
importFrom(utils,capture.output)
useDynLib(bdrc, .registration = TRUE)
14 changes: 13 additions & 1 deletion NEWS.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,16 @@
# bdrc (development version)
# bdrc 2.0.0

* Integrated C++ via Rcpp and RcppArmadillo packages for significant performance improvements.
* Multiple functions rewritten in C++ to speed up the MCMC sampling algorithm and various other tasks.
* The "Deviance" posterior output has been transformed and renamed "Posterior log-likelihood". The Deviance was previously calculated as -2 times the Posterior log-likelihood.
* The plot(tournament_obj, type = "Deviance") figure is now created by evaluating plot(tournament_obj, type = "boxplot").
* The log-likelihood of the models is now computed with the log-transformed discharge observations (normally distributed), rather than with discharge on the real scale (log-normally distributed).
* Pointwise WAIC values (WAIC_i) stored to the model objects.
* Implemented standard error computations for WAIC and Delta_WAIC estimates.
* Applied log-sum-exp trick in WAIC and Bayes factor calculations for numerical stability.
* Added DOI links to references.
* Revised summary() output for tournament objects.
* Updated package vignettes to reflect recent changes and improvements.

# bdrc 1.1.0

Expand Down
119 changes: 119 additions & 0 deletions R/RcppExports.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,119 @@
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

gplm_density_evaluation_unknown_c_cpp <- function(theta, P, h, B, dist, A, y, epsilon, h_min, nugget, n_unique, mu_x, Sig_ab, Z, lambda_c, lambda_sb, lambda_pb, lambda_eta_1, lambda_seta) {
.Call(`_bdrc_gplm_density_evaluation_unknown_c_cpp`, theta, P, h, B, dist, A, y, epsilon, h_min, nugget, n_unique, mu_x, Sig_ab, Z, lambda_c, lambda_sb, lambda_pb, lambda_eta_1, lambda_seta)
}

gplm_density_evaluation_known_c_cpp <- function(theta, P, h, B, dist, A, y, epsilon, nugget, n_unique, mu_x, Sig_ab, Z, lambda_sb, lambda_pb, lambda_eta_1, lambda_seta, c) {
.Call(`_bdrc_gplm_density_evaluation_known_c_cpp`, theta, P, h, B, dist, A, y, epsilon, nugget, n_unique, mu_x, Sig_ab, Z, lambda_sb, lambda_pb, lambda_eta_1, lambda_seta, c)
}

gplm_predict_u_unknown_c_cpp <- function(theta, x, P, B_u, h_unique, h_u, dist_all, h_min, nugget, n_unique, n_u) {
.Call(`_bdrc_gplm_predict_u_unknown_c_cpp`, theta, x, P, B_u, h_unique, h_u, dist_all, h_min, nugget, n_unique, n_u)
}

gplm_predict_u_known_c_cpp <- function(theta, x, P, B_u, h_unique, h_u, dist_all, c, nugget, n_unique, n_u) {
.Call(`_bdrc_gplm_predict_u_known_c_cpp`, theta, x, P, B_u, h_unique, h_u, dist_all, c, nugget, n_unique, n_u)
}

gplm0_density_evaluation_unknown_c_cpp <- function(theta, h, y, A, dist, epsilon, h_min, nugget, n_unique, mu_x, Sig_ab, Z, lambda_c, lambda_se, lambda_sb, lambda_pb) {
.Call(`_bdrc_gplm0_density_evaluation_unknown_c_cpp`, theta, h, y, A, dist, epsilon, h_min, nugget, n_unique, mu_x, Sig_ab, Z, lambda_c, lambda_se, lambda_sb, lambda_pb)
}

gplm0_density_evaluation_known_c_cpp <- function(theta, h, y, A, dist, epsilon, c, nugget, n_unique, mu_x, Sig_ab, Z, lambda_se, lambda_sb, lambda_pb) {
.Call(`_bdrc_gplm0_density_evaluation_known_c_cpp`, theta, h, y, A, dist, epsilon, c, nugget, n_unique, mu_x, Sig_ab, Z, lambda_se, lambda_sb, lambda_pb)
}

gplm0_predict_u_unknown_c_cpp <- function(theta, x, h_unique, h_u, dist_all, h_min, nugget, n_unique, n_u) {
.Call(`_bdrc_gplm0_predict_u_unknown_c_cpp`, theta, x, h_unique, h_u, dist_all, h_min, nugget, n_unique, n_u)
}

gplm0_predict_u_known_c_cpp <- function(theta, x, h_unique, h_u, dist_all, c, nugget, n_unique, n_u) {
.Call(`_bdrc_gplm0_predict_u_known_c_cpp`, theta, x, h_unique, h_u, dist_all, c, nugget, n_unique, n_u)
}

matMult <- function(A, B) {
.Call(`_bdrc_matMult`, A, B)
}

choleskyDecomp <- function(X) {
.Call(`_bdrc_choleskyDecomp`, X)
}

solveArma <- function(A, B) {
.Call(`_bdrc_solveArma`, A, B)
}

matInverse <- function(A) {
.Call(`_bdrc_matInverse`, A)
}

compute_L <- function(X, Sig_x, Sig_eps, nugget) {
.Call(`_bdrc_compute_L`, X, Sig_x, Sig_eps, nugget)
}

compute_w <- function(L, y, X, mu_x) {
.Call(`_bdrc_compute_w`, L, y, X, mu_x)
}

get_MCMC_summary_cpp <- function(X, h) {
.Call(`_bdrc_get_MCMC_summary_cpp`, X, h)
}

calc_variogram_arma <- function(param_mat, i, burnin, nr_iter) {
.Call(`_bdrc_calc_variogram_arma`, param_mat, i, burnin, nr_iter)
}

variogram_chain <- function(T_max, param_mat1, param_mat2, burnin, nr_iter) {
.Call(`_bdrc_variogram_chain`, T_max, param_mat1, param_mat2, burnin, nr_iter)
}

distance_matrix <- function(x) {
.Call(`_bdrc_distance_matrix`, x)
}

create_A_cpp <- function(h) {
.Call(`_bdrc_create_A_cpp`, h)
}

pri <- function(type, args) {
.Call(`_bdrc_pri`, type, args)
}

chain_statistics_cpp <- function(chains) {
.Call(`_bdrc_chain_statistics_cpp`, chains)
}

plm_density_evaluation_unknown_c_cpp <- function(theta, P, h, B, y, epsilon, Sig_x, mu_x, h_min, nugget, lambda_c, lambda_eta_1, lambda_seta) {
.Call(`_bdrc_plm_density_evaluation_unknown_c_cpp`, theta, P, h, B, y, epsilon, Sig_x, mu_x, h_min, nugget, lambda_c, lambda_eta_1, lambda_seta)
}

plm_density_evaluation_known_c_cpp <- function(theta, P, h, B, y, epsilon, Sig_x, mu_x, c, nugget, lambda_eta_1, lambda_seta) {
.Call(`_bdrc_plm_density_evaluation_known_c_cpp`, theta, P, h, B, y, epsilon, Sig_x, mu_x, c, nugget, lambda_eta_1, lambda_seta)
}

plm_predict_u_unknown_c_cpp <- function(theta, x, P, B_u, h_u, h_min, n_u) {
.Call(`_bdrc_plm_predict_u_unknown_c_cpp`, theta, x, P, B_u, h_u, h_min, n_u)
}

plm_predict_u_known_c_cpp <- function(theta, x, P, B_u, h_u, c) {
.Call(`_bdrc_plm_predict_u_known_c_cpp`, theta, x, P, B_u, h_u, c)
}

plm0_density_evaluation_unknown_c_cpp <- function(theta, h, y, epsilon, Sig_ab, mu_x, h_min, nugget, lambda_c, lambda_se) {
.Call(`_bdrc_plm0_density_evaluation_unknown_c_cpp`, theta, h, y, epsilon, Sig_ab, mu_x, h_min, nugget, lambda_c, lambda_se)
}

plm0_density_evaluation_known_c_cpp <- function(theta, h, y, epsilon, Sig_ab, mu_x, c, nugget, lambda_se) {
.Call(`_bdrc_plm0_density_evaluation_known_c_cpp`, theta, h, y, epsilon, Sig_ab, mu_x, c, nugget, lambda_se)
}

plm0_predict_u_unknown_c_cpp <- function(theta, x, h_u, h_min) {
.Call(`_bdrc_plm0_predict_u_unknown_c_cpp`, theta, x, h_u, h_min)
}

plm0_predict_u_known_c_cpp <- function(theta, x, h_u, c) {
.Call(`_bdrc_plm0_predict_u_known_c_cpp`, theta, x, h_u, c)
}

5 changes: 5 additions & 0 deletions R/bdrc-package.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
## usethis namespace: start
#' @importFrom Rcpp sourceCpp
#' @useDynLib bdrc, .registration = TRUE
## usethis namespace: end
NULL
40 changes: 20 additions & 20 deletions R/data.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Swedish Meteorological and Hydrological Institute.
"norn"
Expand All @@ -16,8 +16,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Swedish Meteorological and Hydrological Institute.
"krokfors"
Expand All @@ -28,8 +28,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Swedish Meteorological and Hydrological Institute.
"spanga"
Expand All @@ -40,8 +40,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Swedish Meteorological and Hydrological Institute.
"skogsliden"
Expand All @@ -52,8 +52,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Swedish Meteorological and Hydrological Institute.
"kallstorp"
Expand All @@ -64,8 +64,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Swedish Meteorological and Hydrological Institute.
"melby"
Expand All @@ -76,8 +76,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Icelandic Meteorological Office, Landsvirkjun - the National Power Company of Iceland, and the Icelandic Road and Coastal Administration.
"jokdal"
Expand All @@ -88,8 +88,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Icelandic Meteorological Office, Landsvirkjun - the National Power Company of Iceland, and the Icelandic Road and Coastal Administration.
"jokfjoll"
Expand All @@ -100,8 +100,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Icelandic Meteorological Office, Landsvirkjun - the National Power Company of Iceland, and the Icelandic Road and Coastal Administration.
"nordura"
Expand All @@ -112,8 +112,8 @@
#'
#' @format A data frame with columns:
#' \describe{
#' \item{W}{Measurements of water stage in meters}
#' \item{Q}{Measurements of water discharge in cubic meters per second}
#' \item{\code{W}}{Measurements of water stage in meters}
#' \item{\code{Q}}{Measurements of water discharge in cubic meters per second}
#' }
#' @source Icelandic Meteorological Office, Landsvirkjun - the National Power Company of Iceland, and the Icelandic Road and Coastal Administration.
"skjalf"
Expand Down
Loading

0 comments on commit dd7eb96

Please sign in to comment.