New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/issue 2462 eff sample size pod types #2477

Closed
wants to merge 8 commits into
base: develop
from

Conversation

Projects
None yet
3 participants
@roualdes
Contributor

roualdes commented Feb 27, 2018

Submission Checklist

  • Run unit tests: ./runTests.py src/test/unit
  • [ x] Run cpplint: make cpplint
  • [ x] Declare copyright holder and open-source license: see below

Summary

Add plain old data structure interface to effective_sample_size.

Intended Effect

New interface for Rstan and Pystan.

How to Verify

New public method effective_sample_size_pod() is added so as to allow unit tests. Run

./runTests.py src/test/unit/mcmc/chains_test.cpp

Side Effects

None intended.

Documentation

None, but happy to add. Don't know where.

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):

California State University, Chico

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

Edward A. Roualdes added some commits Feb 27, 2018

Edward A. Roualdes
add effective_sample_size function for std::vector<std::vector<double…
…> > argument; test with effective_sample_size_pod
Edward A. Roualdes
@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Feb 27, 2018

Thanks! I'll review.

Forget that last response, I thought you were asking me to just create the PR.

@bob-carpenter

Thanks for contributing!

I have a bunch of comments, starting with a question about whether we want the interface to really just be plain old data (POD) types or if we want std::vector instead (not a POD).

The rest of the comments are local comments about doc and code structure.

* @return
*/
double
effective_sample_size(const std::vector<std::vector<double> >& samples)

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

We have to make sure with @bgoodri and @ariddell that this will work for both RStan and PyStan. I was proposing a POD solution with just this signature:

double effective_sample_size(const double* draws, int n_chains, int n_draws_per_chain);  

I don't think we can use 2D arrays because we don't know sizes statically.

where the chains were stored one after the other.

This comment has been minimized.

@roualdes

roualdes Feb 28, 2018

Contributor

Right, a POD solution. I backed away from that after thinking that the original issue had some agreement that the std::vector signature would suffice and finding so many stan::math functions based on std::vectors. But, you're right, confirmation would be good. I'll wait until we hear from @bgoodri and/or @ariddell. Happy to rewrite as you all prefer.

* Current implementation takes the minimum number of samples
* across chains as the number of samples per chain.
*
* @param std::vector<std::vector<double> > samples

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

The parameters and return need textual descriptions. They don't need to have the type duplicated---that's already in the signature.

We try to reserve the "sample" for a set of draws and not use "sample" for a single draw. So this should either be singular "sample" or plural "draws". Or it could be chains.

* across all kept samples.
*
* The implementation is close to the effective sample size
* description in BDA3 (p. 286-287). See more details in Stan

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

"BDA3" isn't an adequate citation. Not everyone's familiar with Andrew's cutesy abbreviations for things. I would just point to the Stan manual here.

* description in BDA3 (p. 286-287). See more details in Stan
* reference manual section "Effective Sample Size".
*
* Current implementation takes the minimum number of samples

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

There's no minimums taken here. It all just uses the cross-chain mean and variance estimates to adjust within chain estimates. But I wouldn't even try to say that, I'd just point to the manual. You can include a URL to the online version of the manual:

http://mc-stan.org/users/documentation

std::vector<std::vector<double> > acov(chains);
for (int chain = 0; chain < chains; chain++) {
stan::math::autocovariance(samples[chain], acov[chain]);

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

This probably got copied from elsewhere, but we don't need stan:: if we're in the stan namespace, so this can just be math::autocovariance.

}
double sum_rho_hat_t = std::accumulate(rho_hat_t.begin(),
rho_hat_t.end(),
0.);

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

[Optional] If you're going to have a 0., I'd suggest 0.0 is more readable.

double sum_rho_hat_t = std::accumulate(rho_hat_t.begin(),
rho_hat_t.end(),
0.);
double ess = chains * n_samples;

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

I don't like this pattern of assigning the wrong value, then correcting it. I think this would be clearer as just:

return chains * n_samples / (1 + 2 * sum_rho_hat_t);
rho_hat_t[t + 2] = rho_hat_t[t + 1];
}
}
double sum_rho_hat_t = std::accumulate(rho_hat_t.begin(),

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

Using _t for a non-typedef will be confusing as that's the convention we use for typedefs. You don't need to put the argument on this, so I'd suggest just changing all rho_hat_t to rho_hat, including sum_rho_hat_t.

for (int chain = 0; chain < num_chains(); chain++) {
Eigen::VectorXd Esample = this->samples(chain, index);
samples[chain].resize(Esample.size());
for (int n = 0; n < Esample.size(); n++) {

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

This can be done with an iterator copy.

std::copy(&Esample[0], &Esample[0] + Esample.size(), samples[chain].begin());

But one of the reasons I was suggesting using POD (plain old data) types is that then you don't need this copy, you can just send in the pointers if you use std::vector<double*> as the type. I'm not quite sure what the best thing to do here is. Given everything else going on, the copies shouldn't be too expensive.
Eventually, we'd want to chase all of this through and change the underlying types of member variables so we don't have to copy. Is the original result stored in this calss in a VectorXd?

double effective_sample_size_pod(const int index) const {
std::vector<std::vector<double> > samples(num_chains());
for (int chain = 0; chain < num_chains(); chain++) {
Eigen::VectorXd Esample = this->samples(chain, index);

This comment has been minimized.

@bob-carpenter

bob-carpenter Feb 27, 2018

Contributor

Local variables shouldn't have upper-case names. We reserve that for types in our own code (not that we've followed this everywhere.

@roualdes

This comment has been minimized.

Contributor

roualdes commented Feb 28, 2018

@bob-carpenter, thank you for your detailed comments. It helps a newbie like me get a feel for the style you're after.

Roughly half of these comments could equally apply throughout chains.hpp. Do you want me to clean up the code within the entire chains.hpp file, or just within the portion of this file that I'm working on?

@ariddell

This comment has been minimized.

Member

ariddell commented Feb 28, 2018

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Feb 28, 2018

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Feb 28, 2018

@ariddell I want to write this without having to do any copying. How about a signature like this:

double effective_sample_size(const std::vector<const double*>& chains);

That will require at most a single pointer copy per chain from whatever the underlying data structure is. For both Eigen::Matrix vectors and row vectors and std::vectorvector types, we can use&x[0]to produce thatconst double*` pointer. Or we could write a generic metaprogram.

We could make it even simpler with

double effective_sample_size(const double**& chains);

but that's actually harder to construct.

I don't think we can get away with just double* as that would force a particular memory layout for sequences of chains on the interfaces.

@ariddell

This comment has been minimized.

Member

ariddell commented Feb 28, 2018

@roualdes

This comment has been minimized.

Contributor

roualdes commented Mar 1, 2018

This might be getting above my pay grade, but here it goes anyway. The problem I see with

double effective_sample_size(const std::vector<const double*>& chains);

or anything that doesn't have a single chain's samples stored in a std::vector is the function stan::math::autocovariance which takes as its first argument a std::vector. Hence, for each chain, we'd have to copy from double* to a std::vector. As far as I can tell, this is exactly what we're trying to avoid.

Instead, what do you say to the following signature?

double effective_sample_size(const std::vector<std::vector<double>* >& chains);

Edward A. Roualdes
working prototype for suggested signature for effective_sample_size(s…
…td::vector<std::vector<double>* >& chains)
@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 2, 2018

My proposal is to fix the external services interfaces so they're as easy to use as possible from the interfaces. So this should really be up to @ariddell and @bgoodri.

Then when we settle on the front-end types, we fix the back-end implementations so they don't involve any copying. This is all just a matter of finagling container types. In fact, if we can support C-style pointer arrays, standard vectors, and Eigen vectors using a template and double operator[](int i). It's easy to map a contiguous block of double values to an Eigen type using Eigen::Map---it doesn't copy.

@ariddell

This comment has been minimized.

Member

ariddell commented Mar 2, 2018

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 2, 2018

I wasn't trying to be that ambitious. I just wanted to query @ariddell and @bgoodri to make sure the interface would be usable by RStan and PyStan. Perhaps the best thing to do will be to have a reference to container type.

template <typename C>
real effective_sample_size(const std::vector<const C&>& chains) const;

where C implements

double operator[](size_t i) const;
size_t size() const;

That requirement to implement size() removes the ability to use plain old data types. If we need to support that, then we need a second argument for the size, too. We could even go totally generic and use something like:

template <typename C>
double effective_sample_size(const C& sample, size_t chains, size_t iterations) const;

where we only need C to be such that sample[i][m] returns the double value of draw m from chain i.

None of these are at all hard to implement---just need to decide on one.

Presumably any of the code can be called. I just want to make it convenient for both RStan and PyStan and CmdStan to call without having to copy the entire posterior sample multiple times to do it. CmdStan can use anything. We can also rewrite the functions to be very generic, as in:

@ariddell

This comment has been minimized.

Member

ariddell commented Mar 2, 2018

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 3, 2018

@ariddell

This comment has been minimized.

Member

ariddell commented Mar 3, 2018

Either one is fine. I prefer the first one:

  • double effective_sample_size(const std::vector<const double*>& chains);
  • double effective_sample_size(const double**& chains);
@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 4, 2018

Thanks, @ariddell---that was the answer we need to code this.

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 5, 2018

@roualdes:

Could you rewrite this with this signature:

double effective_sample_size(const std::vector<const double*>& chains);

No pressure, of course---we can get it done one way or another.

@roualdes

This comment has been minimized.

Contributor

roualdes commented Mar 5, 2018

I'm happy to give it a try. After some research over this weekend, I've got two thoughts going. If you get a minute, any feedback is much appreciated.

  1. We need to accept the number of draws as an argument as well, hence

double effective_sample_size(const std::vector<const double*>& chains, int num_draws);

  1. We need to overload stan::math::autocorrelation and stan::math::autocovariance to play nicely with Eigen::Map(&chains[0], num_draws). This should is possible since Eigen::FFT is overloaded for both std::vector and Eigen::VectorXd -- your LingPipe blog post on autocorrelations calculated with FFTs helped me find this.

As a side benefit, this should allow the current working effective_sample_size to avoid its expensive copying from Eigen::VectorXd to std::vector and back.

If I get 2) working, is there a protocol to create a pull request under stan-dev/math that references this issue?

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 5, 2018

@ariddell

This comment has been minimized.

Member

ariddell commented Mar 5, 2018

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 5, 2018

@ariddell

This comment has been minimized.

Member

ariddell commented Mar 5, 2018

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 5, 2018

@roualdes

This comment has been minimized.

Contributor

roualdes commented Mar 7, 2018

I just pushed my latest efforts, alongside a new stan-dev/math pull request. Some notes:

  1. std::vector does not appear to get copied within Eigen's FFT, however I still chose to go with Eigen data structures.

  2. The signature your after works, but I use many Eigen data structures behind it. Hopefully this is OK, as I believe they are faster than hand written for loops around std::vector.

  3. Eigen::Map does show up, but I don't think any copying is done, what with MatrixBase templates and no assignments to Eigen::VectorXd. You'll find these details on the stan-dev/math side. Please correct me if you disagree about the copying.

  4. In the previous effective_sample_size function there is a comment about unequal numbers of draws across the chains. The new implementation assumes equal number of draws amongst the chains.

  5. I at least tried to correct all the suggestions you had in all previous commits.

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 7, 2018

For (2), Eigen data structures aren't any faster than std::vector in and of themselves (though Eigen has the advantage of not constructing initial membrs, which means it's badly behaved in that it can have uninitialized memory, but it's more efficient). Compound operations like matrix multiplication are much more efficiently written in Eigen than naive loops.

For (3), Eigen Map itself only calls the constructor and copies a pointer. Whenever a Map is assigned to an actual matrix type, a copy will happen. This is true in general with expression templates.

For (4), we never did settle on a way to do ragged chains, so assuming everything's the same size is fine for now.

For (5), thanks!

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Mar 8, 2018

@roualdes

This comment has been minimized.

Contributor

roualdes commented Mar 8, 2018

It's worth mentioning, this won't pass its tests until it can use the code from stan-dev/math#790. However, I respect betanalpha's idea (linked above) that, in a broader context, this sort of code should be less piecemeal. So I'm going to patiently wait on this until directed otherwise. Until then, on to some other "good first issue"s. Thanks.

@roualdes

This comment has been minimized.

Contributor

roualdes commented May 28, 2018

Replaced by #2530

@roualdes roualdes closed this May 28, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment