Commit 8ec8c210 authored by Chris Jewell's avatar Chris Jewell
Browse files

Converted lambda_j/i calculations to use tensorA package.

parent de4c0b08
......@@ -25,7 +25,7 @@ Imports:
gplots,
SPIn,
grDevices
RoxygenNote: 5.0.1
RoxygenNote: 6.0.1
KeepSource: TRUE
LazyData: TRUE
Suggests: testthat
......@@ -10,8 +10,8 @@ Acceptance = R6::R6Class("acceptance",
namesTimes,
namesLocations,
namesTypes,
paramsFix) {
if (paramsFix$alpha == FALSE) {
updateSchema) {
if ('alpha' %in% updateSchema) {
self$alpha <- array(
dim = c(nSources,
nTimes,
......@@ -25,7 +25,7 @@ Acceptance = R6::R6Class("acceptance",
)
}
if (paramsFix$r == FALSE) {
if ('r' %in% updateSchema) {
self$r <- array(
dim = c(nTypes,
nSources,
......
......@@ -8,6 +8,10 @@
#####################################################
#' @importFrom SPIn SPIn
ci_chenShao = function(x, alpha) {
assert_that(is.atomic(alpha),
is.numeric(alpha),
alpha >=0,
alpha <= 1)
n <- length(x)
sorted <- sort(x)
......@@ -36,6 +40,10 @@ ci_chenShao = function(x, alpha) {
ci_percentiles <- function(x, alpha) {
assert_that(is.atomic(alpha),
is.numeric(alpha),
alpha >=0,
alpha <= 1)
n <- length(x)
sorted <- sort(x)
......@@ -51,6 +59,10 @@ ci_percentiles <- function(x, alpha) {
ci_SPIn <- function(x, alpha) {
assert_that(is.atomic(alpha),
is.numeric(alpha),
alpha >=0,
alpha <= 1)
region <- tryCatch({
SPIn::SPIn(x, conf = 1 - alpha)$spin
},
......
This diff is collapsed.
......@@ -8,56 +8,75 @@ Posterior_HaldDP = R6::R6Class(
lambda_i = NULL,
lambda_j = NULL,
lambda_j_prop = NULL,
initialize = function(nSources,
nTimes,
nLocations,
nTypes,
n_iter,
iters = 0,
maxIters = 0,
chunkSize = 2000,
initialize = function(n_iter,
namesSources,
namesTimes,
namesLocations,
namesTypes,
namesIters) {
self$alpha <- array(
dim = c(nSources,
nTimes,
nLocations,
n_iter),
dimnames = list(
source = namesSources,
time = namesTimes,
location = namesLocations,
iter = namesIters
)
)
self$q <- array(
dim = c(
nTypes,
n_iter
),
dimnames = list(
type = namesTypes, iter = namesIters
namesTypes) {
self$maxIters = n_iter
self$alpha = to.tensor(
NA,
dim = list(
'source' = namesSources,
'time' = namesTimes,
'location' = namesLocations,
'iter' = 1:n_iter
)
)
self$q = to.tensor(NA,
dim=list(
'type' = namesTypes,
'iter' = 1:n_iter
))
self$s = to.tensor(NA,
dim=list(
'type' = namesTypes,
'iter' = 1:n_iter
))
self$r = to.tensor(NA,
dim=list(
'type' = namesTypes,
'source' = namesSources,
'time' = namesTimes,
'iter' = 1:n_iter
))
},
sample = function(model) {
# Extend the posterior if we need to
self$iters = self$iters + 1
if (self$iters > self$maxIters)
self$extend(self$chunkSize)
self$s <- array(
dim = c(
nTypes,
n_iter
),
dimnames = list(
type = namesTypes, iter = namesIters
)
)
self$r <- array(
dim = c(nTypes,
nSources,
nTimes,
n_iter),
dimnames = list(
type = namesTypes, source = namesSources, time = namesTimes, iter = namesIters
)
)
# Save chain state
for (time in 1:dim(self$alpha)[2]) {
for (location in 1:dim(self$alpha)[3]) {
self$alpha[, time, location, self$iters] <-
model$alphaNodes[[time]][[location]]$getData()
}
for (source in 1:dim(self$alpha)[1]) {
# TODO: change to an lapply?
self$r[, source, time, self$iters] <-
model$rNodes[[time]][[source]]$getData()
}
}
self$q[, self$iters] <-
model$qNodes$getData()
self$s[, self$iters] <- model$qNodes$s
},
extend = function(iters) {
# Extends the posterior file by iters iterations
rnames = list(1:(self$maxIters + iters))
self$alpha = arrayextend(self$alpha, along = 4, size = iters, newdimnames=rnames)
self$q = arrayextend(self$q, along = 2, size = iters, newdimnames=rnames)
self$s = arrayextend(self$s, along = 2, size = iters, newdimnames=rnames)
self$r = arrayextend(self$r, along = 4, size = iters, newdimnames=rnames)
self$lambda_j = arrayextend(self$lambda_j, along = 4, size = iters, newdimnames=rnames)
self$lambda_i = arrayextend(self$lambda_i, along = 4, size = iters, newdimnames=rnames)
self$lambda_j_prop = arrayextend(self$lambda_j_prop, along = 4, size = iters, newdimnames=rnames)
self$maxIters = self$maxIters + iters
},
calc_lambda_i = function(n_iter,
nTimes,
......@@ -68,31 +87,11 @@ Posterior_HaldDP = R6::R6Class(
namesTypes,
namesIters,
k) {
self$lambda_i <- array(
dim = c(nTypes,
nTimes,
nLocations,
n_iter),
dimnames = list(
type = namesTypes,
time = namesTimes,
location = namesLocations,
iter = namesIters
)
)
if (!is.null(self$q)) { # don't try to calculate before any iterations have occured
for (i in 1:n_iter) {
for (times in 1:nTimes) {
for (locations in 1:nLocations) {
self$lambda_i[, times, locations, i] <- self$q[, i] * self$r[, , times, i] %*% (k[, times] * self$alpha[, times, locations, i])
}
}
}
} else {
self$lambda_i <- NULL
}
k = as.tensor(k)
names(k) = c('source','time')
self$lambda_i = self$q * mul.tensor(self$r,i='source',self$alpha * k,j='source',by=c('time','iter'))
self$lambda_i = reorder.tensor(self$lambda_i, c('type','time','location','iter'))
},
calc_lambda_j = function(n_iter,
nSources,
nTimes,
......@@ -102,32 +101,11 @@ Posterior_HaldDP = R6::R6Class(
namesLocations,
namesIters,
k) {
self$lambda_j <- array(
dim = c(nSources,
nTimes,
nLocations,
n_iter),
dimnames = list(
source = namesSources,
time = namesTimes,
location = namesLocations,
iter = namesIters
)
)
if (!is.null(self$q)) { # don't try to calculate before any iterations have occured
for (i in 1:n_iter) {
for (times in 1:nTimes) {
for (locations in 1:nLocations) {
self$lambda_j[, times, locations, i] <- self$alpha[, times, locations, i] * colSums(self$r[, , times, i] * self$q[, i]) * k[, times]
}
}
}
} else {
self$lambda_j <- NULL
}
# tensorA implementation
k = as.tensor(k)
names(k) = c('source','time')
self$lambda_j = self$alpha * mul.tensor(self$r, 'type', self$q, 'type', by='iter') * k
},
calc_lambda_j_prop = function(n_iter,
nSources,
nTimes,
......@@ -137,24 +115,20 @@ Posterior_HaldDP = R6::R6Class(
namesLocations,
namesIters,
k) {
if (is.null(self$lambda_j)) self$calc_lambda_j(n_iter,
nSources,
nTimes,
nLocations,
namesSources,
namesTimes,
namesLocations,
namesIters,
k)
self$lambda_j_prop <- self$lambda_j
for (iter in 1:dim(self$lambda_j)[4]) {
for (times in 1:nTimes) {
for (locations in 1:nLocations) {
self$lambda_j_prop[,times,locations,iter] <-
self$lambda_j[,times,locations,iter] / sum(self$lambda_j[,times,locations,iter])
}
}
}
if (is.null(self$lambda_j))
self$calc_lambda_j(
n_iter,
nSources,
nTimes,
nLocations,
namesSources,
namesTimes,
namesLocations,
namesIters,
k
)
self$lambda_j_prop = to.tensor(apply(self$lambda_j, c('time','location','iter'), function(x) x/sum(x)))
names(self$lambda_j_prop)[1] = 'source'
}
)
)
......@@ -7,15 +7,33 @@
#####################################################
isFiniteInteger = function(a)
isFiniteInteger = function(x)
{
return(is.finite(a) & isTRUE(all.equal(a, as.integer(a))))
is.finite(x) & isTRUE(all.equal(x, as.integer(x)))
}
isFiniteLogical = function(a)
isFiniteLogical = function(x)
{
return(is.finite(a) & is.logical(a))
is.finite(x) & is.logical(x)
}
isFinitePositive = function(a)
isFinitePositive = function(x)
{
return(is.finite(a) & a > 0)
is.finite(x) & x > 0
}
isFinitePositiveInteger = function(x)
{
isFiniteInteger(x) & isFinitePositive(x)
}
arrayextend = function(x, along, size, newdimnames)
{
oldDims = dim(x)
newdim = oldDims
newdim[along] = newdim[along] + size
rnames = dimnames(x)
rnames[along] = newdimnames
newarray = to.tensor(NA, dim=rnames)
args = list(newarray)
args = append(args, lapply(oldDims, function(x) 1:x))
args[[length(args)+1]] = x
do.call('[<-',args)
}
......@@ -68,7 +68,6 @@ priors <-
# Test suite
test_that("HaldDP model construction", {
skip("Skip whole model for now")
set.seed(1)
model <- HaldDP$new(
data = dat,
......@@ -79,8 +78,13 @@ test_that("HaldDP model construction", {
model$fit_params(n_iter = 100,
burn_in = 0,
thin = 1)
sink("/dev/null")
model$update()
sink()
expect_equal_to_reference(model$summary(), "haldDPres.rds")
expect_equal_to_reference(model$summary(), "haldDPres1.rds")
model$update(n_iter=100)
expect_equal_to_reference(model$summary(), "haldDPres2.rds")
expect_equal_to_reference(model$extract('lambda_j'), "haldDPlambdaj.rds")
#model$fit_params(n_iter = 100,
# burn_in = 10,
# thin = 5)
#expect_equal_to_reference(model$summary(), "haldDPres3.rds")
})
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment