diff --git a/.Rprofile b/.Rprofile index 15b877e..cf2babe 100644 --- a/.Rprofile +++ b/.Rprofile @@ -1,7 +1,7 @@ options( help_type="html", repos=c( - CRAN="https://repo.miserver.it.umich.edu/cran/", + CRAN="https://cloud.r-project.org/", kingaa="https://kingaa.github.io" ), useFancyQuotes=FALSE, diff --git a/DESCRIPTION b/DESCRIPTION index 2dc67e9..cae409a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,14 +1,14 @@ Package: circumstance Type: Package Title: pomp parallelized -Version: 0.0.7.0 -Date: 2023-04-14 +Version: 0.0.8.0 +Date: 2023-04-17 Maintainer: Aaron A. King Description: Helper functions for parallelizing pomp computations. Authors@R: c(person(given=c("Aaron","A."),family="King", role=c("aut","cre"),email="kingaa@umich.edu")) URL: https://kingaa.github.io/circumstance/ -Depends: R(>= 4.1.0), pomp(>= 4.6.4) +Depends: R(>= 4.1.0), pomp(>= 5.1) Imports: methods, foreach, grid, grDevices, graphics, utils Suggests: doFuture, doRNG, dplyr, tidyr, ggplot2 Remotes: kingaa/pomp @@ -25,3 +25,4 @@ Collate: 'pfilter.R' 'mif2.R' 'plot_matrix.R' + 'pmcmc.R' diff --git a/NAMESPACE b/NAMESPACE index e516a29..7ef996b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,6 +8,7 @@ export(plot_matrix) exportMethods(continue) exportMethods(mif2) exportMethods(pfilter) +exportMethods(pmcmc) import(methods) importFrom(foreach,"%dopar%") importFrom(foreach,foreach) @@ -35,5 +36,6 @@ importFrom(grid,xaxisGrob) importFrom(grid,yaxisGrob) importFrom(pomp,mif2) importFrom(pomp,pfilter) +importFrom(pomp,pmcmc) importFrom(utils,globalVariables) importFrom(utils,head) diff --git a/R/mif2.R b/R/mif2.R index 0ea55f7..ca7ef53 100644 --- a/R/mif2.R +++ b/R/mif2.R @@ -48,7 +48,9 @@ setMethod( setMethod( "mif2", signature=signature(data = "ANY", starts = "missing"), - definition = pomp::mif2 + definition = function (data, ...) { + pomp::mif2(data,...) + } ) ##' @rdname mif2 diff --git a/R/pfilter.R b/R/pfilter.R index 648ac8f..0560171 100644 --- a/R/pfilter.R +++ b/R/pfilter.R @@ -45,7 +45,9 @@ setMethod( setMethod( "pfilter", signature=signature(data = "ANY", Nrep = "missing"), - definition = pomp::pfilter + definition = function (data, ...) { + pomp::pfilter(data,...) + } ) ##' @rdname pfilter diff --git a/R/pmcmc.R b/R/pmcmc.R new file mode 100644 index 0000000..1d022c4 --- /dev/null +++ b/R/pmcmc.R @@ -0,0 +1,108 @@ +##' Particle Markov chain Monte Carlo in parallel +##' +##' Runs multiple instances of \code{pmcmc} using \code{foreach}. +##' +##' @name pmcmc +##' @rdname pmcmc +##' @importFrom pomp pmcmc +##' @importFrom foreach foreach %dopar% +##' @include pfilter.R mif2.R +##' +##' @param data passed to \code{\link[pomp:pmcmc]{pomp::pmcmc}} +##' @param starts data frame containing parameters at which to begin iterated filtering +##' @param ... all additional arguments are passed to \code{\link[pomp:pmcmc]{pomp::pmcmc}} +##' +##' @seealso \code{\link[pomp:pmcmc]{pomp::pmcmc}}. +##' +NULL + +setGeneric( + "pmcmc", + function (data, starts, ...) + standardGeneric("pmcmc") +) + +##' @rdname pmcmc +##' @export +setMethod( + "pmcmc", + signature=signature(data = "ANY", starts = "data.frame"), + definition = function (data, starts, ...) { + foreach (iter_i=seq_len(nrow(starts)),.combine=c) %dopar% { + pomp::pmcmc(data,params=starts[iter_i,],...) + } -> res + names(res) <- row.names(starts) + attr(res,"doPar") <- get_doPar_info() + res + } +) + +##' @rdname pmcmc +##' @export +setMethod( + "pmcmc", + signature=signature(data = "ANY", starts = "missing"), + definition = function (data, ...) { + pomp::pmcmc(data,...) + } +) + +##' @rdname pmcmc +##' @export +setMethod( + "pmcmc", + signature=signature(data = "pompList", starts = "missing"), + definition = function (data, ...) { + foreach (iter_i=seq_along(data),.combine=c) %dopar% { + pomp::pmcmc(data[[iter_i]],...) + } -> res + names(res) <- names(data) + attr(res,"doPar") <- get_doPar_info() + res + } +) + +##' @rdname pmcmc +##' @export +setMethod( + "pmcmc", + signature=signature(data = "pfilterList", starts = "missing"), + definition = function (data, ...) { + foreach (iter_i=seq_along(data),.combine=c) %dopar% { + pomp::pmcmc(data[[iter_i]],...) + } -> res + names(res) <- names(data) + attr(res,"doPar") <- get_doPar_info() + res + } +) + +##' @rdname pmcmc +##' @export +setMethod( + "pmcmc", + signature=signature(data = "pmcmcList", starts = "missing"), + definition = function (data, ...) { + foreach (iter_i=seq_along(data),.combine=c) %dopar% { + pomp::pmcmc(data[[iter_i]],...) + } -> res + names(res) <- names(data) + attr(res,"doPar") <- get_doPar_info() + res + } +) + +##' @rdname continue +##' @export +setMethod( + "continue", + signature=signature(object = "pmcmcList"), + definition = function (object, ...) { + foreach (iter_i=seq_along(object),.combine=c) %dopar% { + pomp::continue(object[[iter_i]],...) + } -> res + names(res) <- names(object) + attr(res,"doPar") <- get_doPar_info() + res + } +) diff --git a/inst/NEWS b/inst/NEWS index c10503e..e54b334 100644 --- a/inst/NEWS +++ b/inst/NEWS @@ -1,8 +1,14 @@ _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_c_i_r_c_u_m_s_t_a_n_c_e' -_C_h_a_n_g_e_s _i_n '_c_i_r_c_u_m_s_t_a_n_c_e' _v_e_r_s_i_o_n _0._0._6: +_C_h_a_n_g_e_s _i_n '_c_i_r_c_u_m_s_t_a_n_c_e' _v_e_r_s_i_o_n _0._0._8: - • New ‘plot_matrix’ method for making scatterplot matrices. + • New parallel ‘pmcmc’ methods for parallel particle Markov + chain Monte Carlo. + +_C_h_a_n_g_e_s _i_n '_c_i_r_c_u_m_s_t_a_n_c_e' _v_e_r_s_i_o_n _0._0._7: + + • Information on parallel backend is now stored as an attribute + of returned objects. _C_h_a_n_g_e_s _i_n '_c_i_r_c_u_m_s_t_a_n_c_e' _v_e_r_s_i_o_n _0._0._5: diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 5338f36..49e3cf3 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -1,8 +1,13 @@ \name{NEWS} \title{News for package `circumstance'} -\section{Changes in \pkg{circumstance} version 0.0.6}{ +\section{Changes in \pkg{circumstance} version 0.0.8}{ \itemize{ - \item New \code{plot_matrix} method for making scatterplot matrices. + \item New parallel \code{pmcmc} methods for parallel particle Markov chain Monte Carlo. + } +} +\section{Changes in \pkg{circumstance} version 0.0.7}{ + \itemize{ + \item Information on parallel backend is now stored as an attribute of returned objects. } } \section{Changes in \pkg{circumstance} version 0.0.5}{ diff --git a/man/continue.Rd b/man/continue.Rd index 82346d4..f4f60b5 100644 --- a/man/continue.Rd +++ b/man/continue.Rd @@ -1,15 +1,18 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/continue.R, R/mif2.R +% Please edit documentation in R/continue.R, R/mif2.R, R/pmcmc.R \name{continue} \alias{continue} \alias{continue,missing-method} \alias{continue,ANY-method} \alias{continue,mif2List-method} +\alias{continue,pmcmcList-method} \title{Continue an iterative calculation} \usage{ continue(object, ...) \S4method{continue}{mif2List}(object, ...) + +\S4method{continue}{pmcmcList}(object, ...) } \arguments{ \item{object}{the result of an iterative \pkg{pomp} computation} diff --git a/man/pmcmc.Rd b/man/pmcmc.Rd new file mode 100644 index 0000000..6f010e8 --- /dev/null +++ b/man/pmcmc.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pmcmc.R +\name{pmcmc} +\alias{pmcmc} +\alias{pmcmc,ANY,data.frame-method} +\alias{pmcmc,ANY,missing-method} +\alias{pmcmc,pompList,missing-method} +\alias{pmcmc,pfilterList,missing-method} +\alias{pmcmc,pmcmcList,missing-method} +\title{Particle Markov chain Monte Carlo in parallel} +\usage{ +\S4method{pmcmc}{ANY,data.frame}(data, starts, ...) + +\S4method{pmcmc}{ANY,missing}(data, starts, ...) + +\S4method{pmcmc}{pompList,missing}(data, starts, ...) + +\S4method{pmcmc}{pfilterList,missing}(data, starts, ...) + +\S4method{pmcmc}{pmcmcList,missing}(data, starts, ...) +} +\arguments{ +\item{data}{passed to \code{\link[pomp:pmcmc]{pomp::pmcmc}}} + +\item{starts}{data frame containing parameters at which to begin iterated filtering} + +\item{...}{all additional arguments are passed to \code{\link[pomp:pmcmc]{pomp::pmcmc}}} +} +\description{ +Runs multiple instances of \code{pmcmc} using \code{foreach}. +} +\seealso{ +\code{\link[pomp:pmcmc]{pomp::pmcmc}}. +} diff --git a/tests/mif2-02.png b/tests/mif2-02.png index 380f238..46ee43d 100644 Binary files a/tests/mif2-02.png and b/tests/mif2-02.png differ diff --git a/tests/mif2-03.png b/tests/mif2-03.png index c909ac9..965e4e2 100644 Binary files a/tests/mif2-03.png and b/tests/mif2-03.png differ diff --git a/tests/mif2-04.png b/tests/mif2-04.png index 8ef0f36..199de03 100644 Binary files a/tests/mif2-04.png and b/tests/mif2-04.png differ diff --git a/tests/mif2.R b/tests/mif2.R index a7895b1..2a799f1 100644 --- a/tests/mif2.R +++ b/tests/mif2.R @@ -44,6 +44,8 @@ mfs |> mfs |> logLik() +mfs[[1]] |> mif2() + mfs |> continue(Nmif=5) |> traces(pars=c("alpha_1","alpha_3","loglik")) |> diff --git a/tests/mif2.Rout.save b/tests/mif2.Rout.save index 9f9e5bd..b0a582f 100644 --- a/tests/mif2.Rout.save +++ b/tests/mif2.Rout.save @@ -37,7 +37,7 @@ Attaching package: ‘circumstance’ The following objects are masked from ‘package:pomp’: - continue, mif2, pfilter + continue, mif2, pfilter, pmcmc > library(tidyr) > library(dplyr) @@ -104,6 +104,9 @@ Removed 3 rows containing missing values (`geom_line()`). 1 2 3 -480 -483 -477 > +> mfs[[1]] |> mif2() + +> > mfs |> + continue(Nmif=5) |> + traces(pars=c("alpha_1","alpha_3","loglik")) |> @@ -121,12 +124,12 @@ Removed 3 rows containing missing values (`geom_line()`). + as.data.frame() |> + head() .L1 time y1 y2 ess cond.logLik -1 1 1 -4.05 4.781 96.0 -4.19 -2 1 2 1.83 6.273 58.8 -4.69 -3 1 3 -1.32 7.556 36.5 -5.38 -4 1 4 6.64 5.730 30.4 -5.35 -5 1 5 6.58 1.561 121.1 -3.86 -6 1 6 7.54 0.176 80.8 -4.27 +1 1 1 -4.05 4.781 81.1 -4.27 +2 1 2 1.83 6.273 82.8 -4.26 +3 1 3 -1.32 7.556 30.3 -5.31 +4 1 4 6.64 5.730 36.7 -5.25 +5 1 5 6.58 1.561 108.6 -3.99 +6 1 6 7.54 0.176 70.4 -4.53 > > mfs |> + lapply(as_pomp) |> diff --git a/tests/pfilter.Rout.save b/tests/pfilter.Rout.save index a0c698d..6ddc05b 100644 --- a/tests/pfilter.Rout.save +++ b/tests/pfilter.Rout.save @@ -37,7 +37,7 @@ Attaching package: ‘circumstance’ The following objects are masked from ‘package:pomp’: - continue, mif2, pfilter + continue, mif2, pfilter, pmcmc > library(tidyr) > library(dplyr) diff --git a/tests/plot_matrix-1.png b/tests/plot_matrix-1.png index ddff326..dc60be1 100644 Binary files a/tests/plot_matrix-1.png and b/tests/plot_matrix-1.png differ diff --git a/tests/plot_matrix-2.png b/tests/plot_matrix-2.png index 0213796..1973bca 100644 Binary files a/tests/plot_matrix-2.png and b/tests/plot_matrix-2.png differ diff --git a/tests/plot_matrix-3.png b/tests/plot_matrix-3.png index 9306da5..e2591ab 100644 Binary files a/tests/plot_matrix-3.png and b/tests/plot_matrix-3.png differ diff --git a/tests/plot_matrix-4.png b/tests/plot_matrix-4.png index d5f3061..d3cb705 100644 Binary files a/tests/plot_matrix-4.png and b/tests/plot_matrix-4.png differ diff --git a/tests/plot_matrix-5.png b/tests/plot_matrix-5.png index c8309bc..572de57 100644 Binary files a/tests/plot_matrix-5.png and b/tests/plot_matrix-5.png differ diff --git a/tests/plot_matrix.Rout.save b/tests/plot_matrix.Rout.save index b3ca7b0..1461414 100644 --- a/tests/plot_matrix.Rout.save +++ b/tests/plot_matrix.Rout.save @@ -7,6 +7,8 @@ R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. + Natural language support but running in an English locale + R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. @@ -43,7 +45,7 @@ Attaching package: ‘circumstance’ The following objects are masked from ‘package:pomp’: - continue, mif2, pfilter + continue, mif2, pfilter, pmcmc > > png(filename="plot_matrix-%01d.png",res=100) diff --git a/tests/pmcmc-01.png b/tests/pmcmc-01.png new file mode 100644 index 0000000..f2a5d29 Binary files /dev/null and b/tests/pmcmc-01.png differ diff --git a/tests/pmcmc-02.png b/tests/pmcmc-02.png new file mode 100644 index 0000000..f527e22 Binary files /dev/null and b/tests/pmcmc-02.png differ diff --git a/tests/pmcmc-03.png b/tests/pmcmc-03.png new file mode 100644 index 0000000..7062162 Binary files /dev/null and b/tests/pmcmc-03.png differ diff --git a/tests/pmcmc-04.png b/tests/pmcmc-04.png new file mode 100644 index 0000000..26d4117 Binary files /dev/null and b/tests/pmcmc-04.png differ diff --git a/tests/pmcmc.R b/tests/pmcmc.R new file mode 100644 index 0000000..ccd2207 --- /dev/null +++ b/tests/pmcmc.R @@ -0,0 +1,74 @@ +options(digits=3) +png(filename="pmcmc-%02d.png",res=100) + +set.seed(432621217) + +library(circumstance) +library(tidyr) +library(dplyr) +library(ggplot2) +library(doFuture) +library(doRNG) + +registerDoFuture() +registerDoRNG(311289009) + +plan(sequential) + +ou2() -> ou2 + +coef(ou2) |> parmat(3) |> t() |> as.data.frame() -> starts +starts$alpha_3 <- c(-0.1,0.1,0.3) +starts$alpha_1 <- 0.6 + +ou2 |> + window(end=10) |> + pmcmc( + starts=starts, + Np=500,Nmcmc=100, + proposal=mvn_diag_rw(rw.sd=c(alpha_1=0.001,alpha_3=0.001)) + ) -> pms + +pms |> + traces(pars=c("alpha_1","alpha_3","loglik")) |> + plot() + +pms |> + as.data.frame() |> + head() + +pms |> logLik() + +pms[[1]] |> pmcmc() + +pms |> + continue(Nmcmc=10) |> + traces(pars=c("alpha_1","alpha_3","loglik")) |> + plot() + +pms |> + pmcmc(Nmcmc=10) |> + as.data.frame() |> + head() + +pms |> + lapply(as_pomp) |> + concat() |> + pmcmc( + Np=500,Nmcmc=20, + proposal=mvn_diag_rw(rw.sd=c(alpha_1=0.001,alpha_3=0.001)) + ) |> + traces(pars=c("alpha_1","alpha_3","loglik")) |> + plot() + +pms |> + lapply(as,"pfilterd_pomp") |> + concat() |> + pmcmc( + Nmcmc=100, + proposal=mvn_diag_rw(rw.sd=c(alpha_1=0.001)) + ) |> + traces(pars=c("alpha_1","loglik")) |> + plot() + +dev.off() diff --git a/tests/pmcmc.Rout.save b/tests/pmcmc.Rout.save new file mode 100644 index 0000000..ce2859b --- /dev/null +++ b/tests/pmcmc.Rout.save @@ -0,0 +1,143 @@ + +R version 4.2.3 (2023-03-15) -- "Shortstop Beagle" +Copyright (C) 2023 The R Foundation for Statistical Computing +Platform: x86_64-pc-linux-gnu (64-bit) + +R is free software and comes with ABSOLUTELY NO WARRANTY. +You are welcome to redistribute it under certain conditions. +Type 'license()' or 'licence()' for distribution details. + + Natural language support but running in an English locale + +R is a collaborative project with many contributors. +Type 'contributors()' for more information and +'citation()' on how to cite R or R packages in publications. + +Type 'demo()' for some demos, 'help()' for on-line help, or +'help.start()' for an HTML browser interface to help. +Type 'q()' to quit R. + +> options(digits=3) +> png(filename="pmcmc-%02d.png",res=100) +> +> set.seed(432621217) +> +> library(circumstance) +Loading required package: pomp + +Welcome to pomp! + +As of version 4.6, no user-visible pomp function has a name that +includes a dot ('.'). Function names have been changed to replace the +dot with an underscore ('_'). For more information, see the pomp blog: +https://kingaa.github.io/pomp/blog.html. + + +Attaching package: ‘circumstance’ + +The following objects are masked from ‘package:pomp’: + + continue, mif2, pfilter, pmcmc + +> library(tidyr) +> library(dplyr) + +Attaching package: ‘dplyr’ + +The following objects are masked from ‘package:stats’: + + filter, lag + +The following objects are masked from ‘package:base’: + + intersect, setdiff, setequal, union + +> library(ggplot2) +> library(doFuture) +Loading required package: foreach +Loading required package: future +> library(doRNG) +Loading required package: rngtools +> +> registerDoFuture() +> registerDoRNG(311289009) +> +> plan(sequential) +> +> ou2() -> ou2 +> +> coef(ou2) |> parmat(3) |> t() |> as.data.frame() -> starts +> starts$alpha_3 <- c(-0.1,0.1,0.3) +> starts$alpha_1 <- 0.6 +> +> ou2 |> ++ window(end=10) |> ++ pmcmc( ++ starts=starts, ++ Np=500,Nmcmc=100, ++ proposal=mvn_diag_rw(rw.sd=c(alpha_1=0.001,alpha_3=0.001)) ++ ) -> pms +> +> pms |> ++ traces(pars=c("alpha_1","alpha_3","loglik")) |> ++ plot() +> +> pms |> ++ as.data.frame() |> ++ head() + .L1 time y1 y2 ess cond.logLik +1 1 1 -4.05 4.781 116.79 -3.88 +2 1 2 1.83 6.273 38.53 -4.92 +3 1 3 -1.32 7.556 73.81 -4.50 +4 1 4 6.64 5.730 5.77 -6.83 +5 1 5 6.58 1.561 63.59 -4.46 +6 1 6 7.54 0.176 30.84 -5.25 +> +> pms |> logLik() + 1 2 3 +-47.0 -45.6 -45.5 +> +> pms[[1]] |> pmcmc() + +> +> pms |> ++ continue(Nmcmc=10) |> ++ traces(pars=c("alpha_1","alpha_3","loglik")) |> ++ plot() +> +> pms |> ++ pmcmc(Nmcmc=10) |> ++ as.data.frame() |> ++ head() + .L1 time y1 y2 ess cond.logLik +1 1 1 -4.05 4.781 122.47 -3.87 +2 1 2 1.83 6.273 51.46 -4.88 +3 1 3 -1.32 7.556 77.08 -4.29 +4 1 4 6.64 5.730 6.06 -7.22 +5 1 5 6.58 1.561 51.60 -4.95 +6 1 6 7.54 0.176 35.35 -5.19 +> +> pms |> ++ lapply(as_pomp) |> ++ concat() |> ++ pmcmc( ++ Np=500,Nmcmc=20, ++ proposal=mvn_diag_rw(rw.sd=c(alpha_1=0.001,alpha_3=0.001)) ++ ) |> ++ traces(pars=c("alpha_1","alpha_3","loglik")) |> ++ plot() +> +> pms |> ++ lapply(as,"pfilterd_pomp") |> ++ concat() |> ++ pmcmc( ++ Nmcmc=100, ++ proposal=mvn_diag_rw(rw.sd=c(alpha_1=0.001)) ++ ) |> ++ traces(pars=c("alpha_1","loglik")) |> ++ plot() +> +> dev.off() +null device + 1 +>