Re: rstan: extract vs as.array timing (overriding subscripting)

[Ross Boylan]

Redirecting to -dev.

On Wed, 2013-11-27 at 20:57 -0500, Jiqiang Guo wrote:
> I don't think there is a problem here; in fact as.array for stanfit
> calls extract but with parameters other than the default. It is just
> the results returned in different forms that take different amount of
> time.
>
>
> I do not know how or even if something like stanfit[i, j, k] can be
> implemented in R.
The help for Methods has an example of overriding the [] operator:
setMethod("[",
signature(x = "myFrame"),
function (x, i, j, ..., drop = TRUE)
{
S3Part(x) <- callNextMethod()
x@timestamps <- c(Sys.time(), as.POSIXct(x@timestamps))
x
}
)

I'm not sure if anything needs to go in the namespace to support this.

It looks as if the sim$samples[[ic]][[n]] in a stanfit object gives a
vector whose k'th element is the k'th step of chain ic for parameter n.

I guess the default extract is slower because the permutation requires
more work.

I'm not sure if accessing elements through 2 lists is noticeably slower
than accessing array elements.

Ross
>
>
> Jiqiang
>
>
> On Wed, Nov 27, 2013 at 7:11 PM, Ross Boylan <ro...@biostat.ucsf.edu>
> wrote:
> extract is quite slow compared to as.array (13 seconds vs 0.3
> s). I was
> a little surprised, and wonder if there's anything that can be
> done
> about it.
>
> OTOH, extract seems to use less space (see Mb for VCells),
> even though
> it's a more complicated structure. I suppose it must be
> reusing more of
> the stanfit object.
>
> > gc(); system.time(a <- as.array(myfit)); gc()
> used (Mb) gc trigger (Mb) max used (Mb)
> Ncells 1793313 95.8 2908849 155.4 2908849 155.4
> Vcells 294566656 2247.4 439265199 3351.4 374315769 2855.9
> user system elapsed
> 0.320 0.028 0.347
> used (Mb) gc trigger (Mb) max used (Mb)
> Ncells 1793322 95.8 2908849 155.4 2908849 155.4
> Vcells 308591810 2354.4 461308458 3519.6 374315795 2855.9
> > class(a)
> [1] "array"
> > dim(a)
> [1] 1500 3 3116
> > gc(); system.time(a <- extract(myfit)); gc()
> used (Mb) gc trigger (Mb) max used (Mb)
> Ncells 1793316 95.8 2908849 155.4 2908849 155.4
> Vcells 308591801 2354.4 461308458 3519.6 374315795 2855.9
> user system elapsed
> 12.617 0.060 12.718
> used (Mb) gc trigger (Mb) max used (Mb)
> Ncells 1793451 95.8 4953636 264.6 4953636 264.6
> Vcells 312258836 2382.4 484453880 3696.1 374315795 2855.9
> > length(a)
> [1] 18
>
> My interest arises from wanting to do simulations with the
> values of the
> chains, for which I need to extract the coefficients
> repeatedly. If I
> extract the coefficients each time it will be slow; of course
> I can
> cache the results. But it would be nice if the values could
> be accessed
> efficiently, e.g., by stanfit[i, j, k] where i is chain, j is
> step and k
> is coefficient. The grouping of coefficients as in extract
> may be more
> natural for my intended use--I'm not sure what a good
> interface for that
> would be.(*) Maybe an optional argument j to get just one
> step of the
> chain in extract()...
>
> Ross Boylan
>
> (*) I may be doing out of sample predictions, and so some
> individual
> specific effects will not correspond to any coefficients
> estimated, but
> will need to be drawn from the appropriate hiearchical distn.
>
> --
> You received this message because you are subscribed to the
> Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from
> it, send an email to stan-users+...@googlegroups.com.
> For more options, visit
> https://groups.google.com/groups/opt_out.
>
>
> --
> You received this message because you are subscribed to the Google
> Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Ross Boylan]

[Redirected to -dev]

Profiling the code showed that it was spending most of the time in the
unlist method called by sapply in the first line of
fun1 <- function(par) {
sss <- sapply(tidx[[par]], get_kept_samples2, object@sim)
if (is.list(sss)) sss <- do.call(c, sss)
dim(sss) <- c(sum(n_kept), object@sim$dims_oi[[par]])
dimnames(sss) <- list(iterations = NULL)
sss
}
that is part of the extract method for stanfit. This function is
called only if permute is true.

Why that's so slow I don't know. Using lapply instead of sapply runs
fast for that line, but triggers the do.call(c, sss) on the next line,
which is about as slow as the original. However,
fun1 <- function(par) {
sss <- lapply(tidx[[par]], rstan:::get_kept_samples2, object@sim)
if (is.list(sss)) sss <- do.call(cbind, sss)
dim(sss) <- c(sum(n_kept), object@sim$dims_oi[[par]])
dimnames(sss) <- list(iterations = NULL)
sss
}

Runs much faster (since I defined the test function outside of the
package scope I needed to explicitly qualify calls to internal package
functions).

e1 is the original extract; e2 is the new version.
> system.time(t1 <- e1(myfit))
user system elapsed
11.420 0.004 11.463
> system.time(t2 <- e2(myfit))
user system elapsed
0.396 0.008 0.405
!> identical(t1, t2)
[1] TRUE


I think the lines
if (is.list(sss)) sss <- do.call(cbind, sss)
dim(sss) <- c(sum(n_kept), object@sim$dims_oi[[par]])
could be converted to
do.call(cbind, sss)
but perhaps I'm missing some subtlety. As I said, they look as if they
might be a holdover from an earlier implementation.

Ross

On Wed, 2013-11-27 at 20:57 -0500, Jiqiang Guo wrote:
> I don't think there is a problem here; in fact as.array for stanfit
> calls extract but with parameters other than the default. It is just
> the results returned in different forms that take different amount of
> time.
>
>
> I do not know how or even if something like stanfit[i, j, k] can be implemented in R.
>
>

[Ross Boylan]

On Tue, 2013-12-03 at 14:01 -0800, Ross Boylan wrote:
> Profiling the code showed that it was spending most of the time in the
> unlist method called by sapply in the first line of
> fun1 <- function(par) {
> sss <- sapply(tidx[[par]], get_kept_samples2,
> object@sim)
> if (is.list(sss)) sss <- do.call(c, sss)
> dim(sss) <- c(sum(n_kept), object@sim$dims_oi[[par]])
> dimnames(sss) <- list(iterations = NULL)
> sss
> }
> that is part of the extract method for stanfit. This function is
> called only if permute is true.

The sense in which unlist is called by sapply may not be obvious (it
wasn't to me)--it has nothing to do with the get_kept_samples2 function.
Rather, sapply is
> sapply
function (X, FUN, ..., simplify = TRUE, USE.NAMES = TRUE)
{
FUN <- match.fun(FUN)
answer <- lapply(X = X, FUN = FUN, ...)
if (USE.NAMES && is.character(X) && is.null(names(answer)))
names(answer) <- X
if (!identical(simplify, FALSE) && length(answer))
simplify2array(answer, higher = (simplify == "array"))
else answer
}
and simplify2array ultimately calls unlist.

Here's a sample call stack from Rprof:
"unlist" "as.vector" "simplify2array" "sapply" "FUN" "lapply" ".local" "extract" "standardGeneric" "extract"
(outermost call on the right, apparently, and "FUN" is the function fun1 in the source).

If anyone is curious, here is the profile of the original function:
> Rprof()
> e <- extract(myfit)
> Rprof(NULL)
> summaryRprof()
$by.self
self.time self.pct total.time total.pct
"unlist" 8.30 92.22 8.30 92.22
"FUN" 0.44 4.89 8.90 98.89
"lapply" 0.08 0.89 8.98 99.78
":" 0.04 0.44 0.04 0.44
"array" 0.04 0.44 0.04 0.44
"as.vector" 0.02 0.22 8.32 92.44
"<Anonymous>" 0.02 0.22 0.04 0.44
"do.call" 0.02 0.22 0.04 0.44
"dimnames<-" 0.02 0.22 0.02 0.22
"sub" 0.02 0.22 0.02 0.22

$by.total
total.time total.pct self.time self.pct
"lapply" 8.98 99.78 0.08 0.89
"extract" 8.98 99.78 0.00 0.00
".local" 8.98 99.78 0.00 0.00
"standardGeneric" 8.98 99.78 0.00 0.00
"FUN" 8.90 98.89 0.44 4.89
"sapply" 8.86 98.44 0.00 0.00
"simplify2array" 8.36 92.89 0.00 0.00
"as.vector" 8.32 92.44 0.02 0.22
"unlist" 8.30 92.22 8.30 92.22
":" 0.04 0.44 0.04 0.44
"array" 0.04 0.44 0.04 0.44
"<Anonymous>" 0.04 0.44 0.02 0.22
"do.call" 0.04 0.44 0.02 0.22
"dimnames<-" 0.02 0.22 0.02 0.22
"sub" 0.02 0.22 0.02 0.22
"doTryCatch" 0.02 0.22 0.00 0.00
".ess_funargs" 0.02 0.22 0.00 0.00
"gsub" 0.02 0.22 0.00 0.00
"methods" 0.02 0.22 0.00 0.00
"tryCatch" 0.02 0.22 0.00 0.00
"tryCatchList" 0.02 0.22 0.00 0.00
"tryCatchOne" 0.02 0.22 0.00 0.00

Oddly, the time reported by Rprof was consistently less than reported by
system.time running the unprofiled code: c 9 vs 12 seconds.
Ross

[Jiqiang Guo]

Good finding. Thanks a lot. So I will change those two lines to

sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
sss <- do.call(cbind, sss)

-- Jiqiang

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Jiqiang Guo]

Thanks.  I think it might not be good to do that for 

* it need do the boundary check, which makes it inefficient

* the thinning makes the index for iteration ambiguous 

--

Jiqiang 

You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Ross Boylan]

On Tue, 2013-12-03 at 19:17 -0500, Jiqiang Guo wrote:
> Thanks. I think it might not be good to do that for
> * it need do the boundary check, which makes it inefficient
> * the thinning makes the index for iteration ambiguous

More generally, it's unclear exactly what interface would be of general
interest. I originally was thinking

stanfit[i, j, k] where i is chain, j is step and k
is coefficient.

As you point out, one might want j to be before or after thinning.
k could be an index (in the total index sense used internally), a string
like the values of pars, or booleans. The total index seems an
unfriendly interface, though perhaps if some of the related functions
were exposed it would make sense.

One might want to overlook the chains and use
stanfit[n, k] where n is the n'th sample (combining all chains), and k
indexes a parameter.

It might need some caching (e.g., of the total indices for a named
parameter) to be effficient. I'm not sure it needs bounds checking; I
think it will just throw an error if one attempts an out of bounds
access.

There are also several choices about what gets returned.

All of which is to agree that the time is not ripe for implementing [].
Now that extract is faster I feel less need to work around it.

Ross

> +unsub...@googlegroups.com.

[Ross Boylan]

On Tue, Dec 03, 2013 at 07:14:48PM -0500, Jiqiang Guo wrote:
> Good finding. Thanks a lot. So I will change those two lines to
>
> sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
> sss <- do.call(cbind, sss)
>
> -- Jiqiang
>

Apparently the omitted lines do something, since without them I've
changed the interface for 2 dimensional coefficients, e.g.,
multinomial models.

I've been running with
fun1 <- function(par) {

sss <- lapply(tidx[[par]], get_kept_samples2,
object@sim)
sss <- do.call(cbind, sss)

dimnames(sss) <- list(iterations = NULL)
sss
}

I'll have to think about how to fix it. Simply restoring the
statement setting the dimensions (and probably names, since they
vanished too) may result in some transpositions, and there may be a
faster/better way to combine when tidx[[par]] is a matrix.

Not quite ready for prime-time!

Ross

[Ross Boylan]

On Thu, Dec 05, 2013 at 07:46:36PM -0800, Ross Boylan wrote:
> On Tue, Dec 03, 2013 at 07:14:48PM -0500, Jiqiang Guo wrote:
> > Good finding. Thanks a lot. So I will change those two lines to
> >
> > sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
> > sss <- do.call(cbind, sss)
> >
> > -- Jiqiang
> >
>
> Apparently the omitted lines do something, since without them I've
> changed the interface for 2 dimensional coefficients, e.g.,
> multinomial models.
>
> I've been running with
> fun1 <- function(par) {
> sss <- lapply(tidx[[par]], get_kept_samples2,
> object@sim)
> sss <- do.call(cbind, sss)
> dimnames(sss) <- list(iterations = NULL)
> sss
> }
>
> I'll have to think about how to fix it. Simply restoring the
> statement setting the dimensions (and probably names, since they
> vanished too) may result in some transpositions, and there may be a
> faster/better way to combine when tidx[[par]] is a matrix.
>
> Not quite ready for prime-time!
>
> Ross

Adding back the deleted dim-setting line seems to work:

fun1 <- function(par) {
sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
sss <- do.call(cbind, sss)

# essential for 2 dimensional parameters

[Ross Boylan]

On 12/6/2013 8:22 AM, Ross Boylan wrote:
> On Thu, Dec 05, 2013 at 07:46:36PM -0800, Ross Boylan wrote:
>> On Tue, Dec 03, 2013 at 07:14:48PM -0500, Jiqiang Guo wrote:
>>> Good finding. Thanks a lot. So I will change those two lines to
>>>
>>> sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
>>> sss <- do.call(cbind, sss)
>>>
>>> -- Jiqiang
>>>
>>

> Adding back the deleted dim-setting line seems to work:
> fun1 <- function(par) {
> sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
> sss <- do.call(cbind, sss)
> # essential for 2 dimensional parameters
> dim(sss) <- c(sum(n_kept), object@sim$dims_oi[[par]])
> dimnames(sss) <- list(iterations = NULL)
> sss
> }

BTW, on a quick test the code is a bit over 10x faster than the current
production code for 2-dimensional coefficients.

[Jiqiang Guo]

I am a little bit confused.  Anyway, the is the code in current develop branch,

https://github.com/stan-dev/rstan/blob/develop/rstan/rstan/R/stanfit-class.R#L384

I have never left out the code for setting the dimensions. 

--

Jiqiang 

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Ross Boylan]

On 12/6/2013 10:57 AM, Jiqiang Guo wrote:

I am a little bit confused.  Anyway, the is the code in current develop branch,

https://github.com/stan-dev/rstan/blob/develop/rstan/rstan/R/stanfit-class.R#L384

I have never left out the code for setting the dimensions.

I did that in my local copy.  I wasn't sure what you did from your earlier note.
So you have no problem to fix.

Ross

--

Jiqiang 

On Fri, Dec 6, 2013 at 12:53 PM, Ross Boylan <ro...@biostat.ucsf.edu> wrote:

On 12/6/2013 8:22 AM, Ross Boylan wrote:

On Thu, Dec 05, 2013 at 07:46:36PM -0800, Ross Boylan wrote:

On Tue, Dec 03, 2013 at 07:14:48PM -0500, Jiqiang Guo wrote:

Good finding. Thanks a lot. So I will change those two lines to

sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
sss <- do.call(cbind, sss)

-- Jiqiang

Adding back the deleted dim-setting line seems to work:
             fun1 <- function(par) {
               sss <- lapply(tidx[[par]], get_kept_samples2, object@sim)
               sss <- do.call(cbind, sss)
               # essential for 2 dimensional parameters
               dim(sss) <- c(sum(n_kept), object@sim$dims_oi[[par]])
               dimnames(sss) <- list(iterations = NULL)
               sss
             }

BTW, on a quick test the code is a bit over 10x faster than the current production code for 2-dimensional coefficients.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.