help with using stan::mcmc::chains

[Jiqiang Guo]

I need some help with using stan::mcmc::chains.  My current C++ skills 

could not solve it though I am still looking for solutions from the books. 

I would like to have a member of class stan::mcmc::chains in a class, which would be

exposed to R by Rcpp's module. I could initialize the chains object in the class' 

constructor initialization list. That is, suppose my new class is stan_fit, I have 

template <Model m>
class stan_fit {
  stan::mcmc::chains c_;
  stan_fit(SEXP data, SEXP args) : c_(...) 

  {}
};

However, in order to construct a chains object, I need the parameter names and dimensions,
which have to be get by call functions of the model. Then I guess I cannot use 

the initialization list.  Nor could I use assign operation to reassign the chains object to a new 

chains object inside the ctor as follows (I guess because there are const members of class chains). 

std::vector<std::string> pnames;

std::vector<std::vector<size_t> > pdimss;

model.get_param_names(pnames);

model.get_dims(pdimss);

c_ = stan::mcmc::chains(num_chains, pnames, pdimms); // seems to be problematic

Thanks,

Jiqiang 

[Bob Carpenter]

There are some design questions down below about
the model object that others might be interested in.

On 5/26/12 4:41 PM, Jiqiang Guo wrote:
> I need some help with using stan::mcmc::chains. My current C++ skills
> could not solve it though I am still looking for solutions from the books.

Your best bet for this kind of thing is the
Myers book, but it's hardly comprehensive. And
it's a hard sort of thing to look up.

> I would like to have a member of class stan::mcmc::chains in a class, which would be
> exposed to R by Rcpp's module. I could initialize the chains object in the class'
> constructor initialization list. That is, suppose my new class is stan_fit, I have
>
> template <Model m>
> class stan_fit {
> stan::mcmc::chains c_;
> stan_fit(SEXP data, SEXP args) : c_(...)
> {}
> };
>
> However, in order to construct a chains object, I need the parameter names and dimensions,
> which have to be get by call functions of the model. Then I guess I cannot use
> the initialization list. Nor could I use assign operation to reassign the chains object to a new
> chains object inside the ctor as follows (I guess because there are const members of class chains).
>
> std::vector<std::string> pnames;
> std::vector<std::vector<size_t> > pdimss;
> model.get_param_names(pnames);
> model.get_dims(pdimss);
> c_ = stan::mcmc::chains(num_chains, pnames, pdimms); // seems to be problematic

I have a tendency to make everything I can immutable (const) and
set in the constructor rather than with setters/getters. The
main advantage is that it's easier to reason about the state of
a class -- it's (a) thread safe for access to the immutables upon
construction, and (b) the immutables are always well-defined upon
construction, meaning less state checking for error handling (like
what get_names() should do before the names have been set).

It can lead to some clunky constructors, though.

First, here's how to code around the problem.
You can write a "helper" function to get the name vector:

template <class Model>
std::vector<std::string> get_names(const Model& m) {
std::vector<std::string> names;
m.get_param_names(names);
return names; // copy for return
}

and similarly a get_dimss() for the dimensions. Then call them
in the constructor:

template <class Model>
class stan_fit {

private:
const Model& _model; // store externally, hold ref
stan::mcmc::chains _chains; // store locally, hold obj

public:
stan_fit(const Model& model, size_t num_chains)
: _model(model),
_chains(num_chains,
get_names(model), // copy for model
get_dimss(model)) {
}

...

};

I think it might be easier if we make a few changes:

1. Have the models extend an abstract base class that itself
extends the prob-gradient class with all the methods you'll
need. Then you don't need to template stan_fit on a Model
concept.

2. Have the get_param_names() function return a vector. Then
you won't need the helper function. I can't see many operations
where this extra copy's going to be a bottleneck. Or,
we can implement the helper function in the abstract base class.

- Bob