proposal draft for RStan interface

[Bob Carpenter]

Ahead of working on the Python/R under-the-hood refactor when
Allen gets here, I tried to go through the current R interface and
whittle things down (and expand a bit) to where I think it should go.
I'd think that roughly the same interface would work in Python.

I hope this is close enough to idiomatic for R that it's not
too embarrassing --- but please send feedback (up to and including
telling me I should just stay out of R!).

The biggest change is simplifying what functions do and where some
of the action is, as well as splitting some of the functions out (like
stan() is now stan_nuts, stan_hmc, and stan_diagnose).

I don't think this loses any functionality other than the big fit
object, which I don't think we should have. And it's assuming some
functionality we don't yet have, like being able to input a mass matrix.

I can move this over to a wiki (on RStan or Stan's repo?) if people
think that'll be easier and the response isn't to just get out of the
way of RStan.

- Bob


1. Compilation
------------------------------------------------------------
The compilation function compiles a Stan program to a dynamic shared
object (DSO) and optionally saves it.

compiled_program
rtan::stan_compile( # DEFAULTS
source = (<string> || <path>) # no default
name = <string>, # empty
boost = <path>, # use BH
eigen = <path>, # use RcppEigen
opt = ("fast" || "debug"), # "fast"
save_dso = <bool>, # FALSE
debug_code = <bool> # TRUE
);

* Is there any reason to keep the name option?
* How long-lived is the "save_dso" option?
- Where is it saved?
* If we allow optimize settings, is there a way to have it scope
over just the call and not affect the R global build environment?
* Is there a way we can have make-like behavior here where if there's
a saved version and the model file or code have not changed, don't
do anything?
* Do we need the other existing compiler options ("presentation1",
"presentation2", and "small")?
* Should the save_dso provide a path rather than a boolean?
* Is it possible to save the DSO across sessions?
* Should the dave_dso defalt be true and if so, use temp?
* Do we ever need NDEBUG=false? That is, shouldn't we just always
have -DNDEBUG set to turn off debugging? (The double-negative

compiled_program
--------------------
This should only include the compiled Stan program and methods
to do variable conversion, extract names, compute log densities, etc.

Slots:
name : <string>
stan_code : <string>
cpp_code : <string>
dso : <S4 cxxdso> (from inline)
cxx_flags: <string>
Methods:
double
log_prob(
unconstrained_params = <vector>,
jacobian = <bool>);

list(value=<double>, grad=<vector>)
log_prob_grad(
unconstrained_params = <vector>,
jacobian = <bool>);

return$value : <double>
return$grad : <vector>

list(
constrain_params(
unconstrained_params = <vector>);

vector
constrain_params_vec(
unconstrained_params = <vector>);

vector
unconstrain_params(
params = (<list> || vector));

* Seems there should be a better way to organize this, but I
can't see it. And maybe some shorter names


2. NUTS
------------------------------------------------------------

mcmc_fit
rstan::stan_nuts(
program = <compiled_program>, # no default
data = (<list> || <environment>), # top-level environment by default
max_treedepth = <int>, # 10
chain_ids = <ints>, # 1:4
num_warmup = <int>, # 10
num_samples = <int>, # 10
init=(<double> # 2.0
|| <list>
|| <list of lists>
|| <function>)
stepsize=<double>, # 1.0
metric = ("unit_e", # "diag_e"
|| "diag_e"
|| "dense_e"
|| <vector>
|| <matrix>),
adapt_stepsize = <bool>, # TRUE
adapt_mass_matrix = <bool>, # TRUE
adapt = list(gamma=<double>,
kappa=<double>,
t0=<double>,
init_buffer=<int>,
term_buffer=<int>,
adapt_window=<int>) # ...
params_to_save = (NA
|| <vector(char)>) # NA
save_warmup = <bool>, # TRUE
thin = <int>, # 1
sample_file = <path>, # NA
diagnostic_file = <path>, # NA
refresh = <int>, # 1
verbose = ("quiet"
|| "standard"
|| "verbose"), # standard
seed = <int> # random from R
)

* Return value just includes the result in such a way that it
can be extracted as before.
* Return value won't include the compiled program (or source).
* Should return value include timing information?
* There's no append_samples = <bool>
* There's no option to compile the model from within the stan_nuts
call.
* data defaults to top-level environment with lazy read (not how
it's done now); list is list of values as before
* Instead of "chains" and a start ID, there's just a "chain_ids"
argument supplying both the chains and the ID, defaulting to the usual
* Default num_warmup and num_samples very small by default (but shouldn't
be so small as to trigger scary messages about adaptation)
* Metric is specified by name of type or by value; the names default
to unit values for init
* Those adapt parameters are all for adapt_stepsize, right?
* Are there adaptation parameters for mass_matrix, like regularization?
* params_to_save defaults to NA (or all), or a positive list of whole
parameter names.
* sample_file gets a copy of the output written to it, but return
object still created --- should this just move so that the return
object has a save on it? Or should setting this cause the return of
this function to be NA?

mcmc_fit
--------------------


3. EHMC
------------------------------------------------------------
Plain old HMC looks just like NUTS, but with max_treedepth replaced
by integration time:

mcmc_fit
rstan::stan_hmc(...
int_time = <double>
...
)


4. Optimize
------------------------------------------------------------

optimization_fit
stan_lbfgs(
program = <compiled_program>,
data = (<list> || <environment>),
init = (<double>
|| <list>
|| <list of lists>
|| <function>),
optimization_file = <file>,
verbose = ("quiet" || "verbose"),
hessian = <bool>, # FALSE
max_iterations = <int>, # NA
save_iterations = <bool>, # FALSE
refresh = <int>, # 5
init_alpha = <double>, # ???
tol_obj = <double>,
tol_grad = <double>,
tol_param = <double>,
tol_obj_rel = <double>,
tol_grad_rel = <double>,
tol_param_rel = <double>,
seed = <int>
)


* Should verbose be a boolean?
* What does current as_vector do? Current doc says "flag indicating whether a
vector is used to for the point estimate found. A list can be
used instead by specifying it to be FALSE"
* What's the behavior of tolerances if more than one's specified? I'd
think it should be read as AND or as OR.
* No option to specify parameter return as vector --- use a compiled
model's conversion.

optimization_fit
--------------------
$par : parameters (as list)
$value : lp__


5. Diagnose
------------------------------------------------------------

diagnostics
rstan::stan_diagnose(
program = <compiled_program>,
data = (<list> || <environment>),
init = (<double>
|| <list>
|| <list of lists>
|| <function>),
epsilon = <double>,
error_threshold = <double>,
seed = <int>,
chain_ids = <int>)

* arguments just a subset of those for NUTS/HMC
* returns diagnostics for each of the chain_ids
* not sure what goes in diagnostics as long as it's
easy to print

diagnostics
--------------------
?? not sure what should go in here --- it needs to be printable

==============================================================================

[Allen B. Riddell]

Hi Bob,

This is great. It looks like you incorporated all the ideas from the
wiki (https://github.com/stan-dev/rstan/wiki/Interfaces-3.0-Spec).
Perhaps the RStan wiki is the right place for this discussion, since I
think we'll want to be making minor modifications as we go along and
inviting more comments.

The essentials for me are:

1. Consistent variable names across Stan/PyStan/RStan (e.g., `warmup` in
RStan/PyStan needs to become `num_warmup`).
2. Intermediate classes/objects that are more self-contained.

Both of these items are covered by your proposal/email, so my comments
are mostly concerning details -- I'll put them inline below. I'd also be
happy to make the edits on the wiki.

On 05/21, Bob Carpenter wrote:
> The biggest change is simplifying what functions do and where some
> of the action is, as well as splitting some of the functions out (like
> stan() is now stan_nuts, stan_hmc, and stan_diagnose).
> I don't think this loses any functionality other than the big fit
> object, which I don't think we should have. And it's assuming some
> functionality we don't yet have, like being able to input a mass matrix.
>

Sounds good.

>
> 1. Compilation
> ------------------------------------------------------------
> The compilation function compiles a Stan program to a dynamic shared
> object (DSO) and optionally saves it.
>
> compiled_program
> rtan::stan_compile( # DEFAULTS
> source = (<string> || <path>) # no default
> name = <string>, # empty
> boost = <path>, # use BH
> eigen = <path>, # use RcppEigen
> opt = ("fast" || "debug"), # "fast"
> save_dso = <bool>, # FALSE
> debug_code = <bool> # TRUE
> );
>
> * Is there any reason to keep the name option?
> * How long-lived is the "save_dso" option?
> - Where is it saved?

I think we could get rid of both name and save_dso. Whatever it is that
is returned (`compiled_program`) can be serialized and saved if one
wants to save it. R and Python both make this easy. Disk space is not a
concern. If we need a unique identifier for the model we can take the
hash of the generated C++ source code.

> * If we allow optimize settings, is there a way to have it scope
> over just the call and not affect the R global build environment?
> * Is there a way we can have make-like behavior here where if there's
> a saved version and the model file or code have not changed, don't
> do anything?

If a compiled program's identifier is the hash of the generated C++ code
we have any easy identity check!

> * Do we need the other existing compiler options ("presentation1",
> "presentation2", and "small")?
> * Should the save_dso provide a path rather than a boolean?
> * Is it possible to save the DSO across sessions?
> * Should the dave_dso defalt be true and if so, use temp?
> * Do we ever need NDEBUG=false? That is, shouldn't we just always
> have -DNDEBUG set to turn off debugging? (The double-negative
>

It would be great to clean up the flags.

re: the function name, perhaps just rstan::compile/pystan.compile?

LGTM

Yes! This is essential to making serialization easy. And given how long
it takes to compile programs being more nimble with serializing and
caching models seems like it might be important in the future.

> * Should return value include timing information?
> * There's no append_samples = <bool>
> * There's no option to compile the model from within the stan_nuts
> call.

Yes.

> * data defaults to top-level environment with lazy read (not how
> it's done now); list is list of values as before
> * Instead of "chains" and a start ID, there's just a "chain_ids"
> argument supplying both the chains and the ID, defaulting to the usual
> * Default num_warmup and num_samples very small by default (but shouldn't
> be so small as to trigger scary messages about adaptation)
> * Metric is specified by name of type or by value; the names default
> to unit values for init
> * Those adapt parameters are all for adapt_stepsize, right?
> * Are there adaptation parameters for mass_matrix, like regularization?
> * params_to_save defaults to NA (or all), or a positive list of whole
> parameter names.
> * sample_file gets a copy of the output written to it, but return
> object still created --- should this just move so that the return
> object has a save on it? Or should setting this cause the return of
> this function to be NA?

re: sample_file, for simplicity I think we should get rid of it here and add
a `to_csv` method or similar to the mcmc_fit object. If people are concerned
about memory they should be using cmdstan.

---

LGTM in general. It's an implementation detail but I would prefer to
have this as a method of compiled_program (in R and Python). Perhaps
`sample` could be an alias to this nuts method? Could we call it `nuts`
instead of `stan_nuts`?

> mcmc_fit
> --------------------
>

right now the method is `extract` -- we should think about what this
could be. I wonder if we could do some black magic and add some lazy
extraction methods. For example, why not have mcmc_fit['mu'] (Python) or
mcmc_fit$mu extract the samples?

>
> 3. EHMC
> ------------------------------------------------------------
> Plain old HMC looks just like NUTS, but with max_treedepth replaced
> by integration time:
>
> mcmc_fit
> rstan::stan_hmc(...
> int_time = <double>
> ...
> )

Sounds good.

as_vector was a feature that was added later in the life of Stan 2.0. I
think we should pick a default. It's basically letting the user decide
between a return value that is a param-value map (list/dictionary) or a
vector (with names). I'd vote for the vector being the default,
especially in R where this is trivial. If people want the `lp__` at that
point they should make an extra call to `log_prob`.

> * What's the behavior of tolerances if more than one's specified? I'd
> think it should be read as AND or as OR.
> * No option to specify parameter return as vector --- use a compiled
> model's conversion.
>

general comments: again, can this be a method of compiled_program? Can optimize be an
alias?

Same comment about optimization_file as for sample_file above.

> optimization_fit
> --------------------
> $par : parameters (as list)
> $value : lp__
>

I think it should just be parameters (with names) as a vector. If people
want the lp__ value they can get it as an extra step by calling
`log_prob`.

>
> 5. Diagnose
> ------------------------------------------------------------
>
> diagnostics
> rstan::stan_diagnose(
> program = <compiled_program>,
> data = (<list> || <environment>),
> init = (<double>
> || <list>
> || <list of lists>
> || <function>),
> epsilon = <double>,
> error_threshold = <double>,
> seed = <int>,
> chain_ids = <int>)
>
> * arguments just a subset of those for NUTS/HMC
> * returns diagnostics for each of the chain_ids
> * not sure what goes in diagnostics as long as it's
> easy to print
>
> diagnostics
> --------------------
> ?? not sure what should go in here --- it needs to be printable
>

Can this be a method and be called rstan::diagnose?

Thanks,

Allen

[Michael Betancourt]

So my big problem with this flat structure is a conceptual one.
Right now there’s a glut of terrible algorithm names out there
because people don’t understand the foundations of HMC.
Consequently "Hamiltonian Monte Carlo” is often relegated
to what is more properly called Metropolis-Corrected Euclidean
Hamiltonian Monte Carlo with a static integration time, and
every tweak yields a new name. This flurry of names obfuscates
the real underlying mathematical framework — changing the
metric doesn’t really lead to a new algorithm, just a different
implementation. Why, for example, would RHMC receive
its own command but EHMC with a diagonal metric vs EHMC
with a dense metric would be different options within the same
command? Moreover, a single command “stan_hmc” that
just calls one of this configurations makes everything so much
worse.

At the very least we need to be more careful with the names.
I know we all love NUTS, but it’s really just an implementation
of a more general idea of dynamic integration time (really the
no-U-turn part of it is the least novel — the really important
part of the algorithm is binary tree to preserve reversibility).
What we have are two main different configuration options:

Metric: Identity, diagonal Euclidean, dense Euclidean (and in the future SoftAbs Riemannian)
Integration Time: Static vs Exhaustive (i.e. integrate longer on bigger level sets, NUTS is an example)

Why is integration time a better split than metric?

Then there’s the fact that we can sample or optimize, and
the current names make no distinction between the two
If people are that against hierarchical arguments (note
that we have yet to have any complaint from users on
the list — the CmdStan user base is small, but it’s not
negligible) then we at least have to maintain the hierarchy
is the command names: stan_mcmc_*, stan_opt_*, stan_vb_*, etc.

Finally, maintaining this orgy of commands is going to be a
pain in the ass. What to change one argument that all of the
commands share? Go ahead and tweak every single command
config.

> --
> 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/d/optout.

[Allen B. Riddell]

These are good points. And they relate to the design of the "config
object" as well as the organization of a family of functions.

> Then there’s the fact that we can sample or optimize, and
> the current names make no distinction between the two
> If people are that against hierarchical arguments (note
> that we have yet to have any complaint from users on
> the list — the CmdStan user base is small, but it’s not
> negligible) then we at least have to maintain the hierarchy
> is the command names: stan_mcmc_*, stan_opt_*, stan_vb_*, etc.
>
> Finally, maintaining this orgy of commands is going to be a
> pain in the ass. What to change one argument that all of the
> commands share? Go ahead and tweak every single command
> config.
>

I agree that maintaining this multitude of commands is going to be a
pain in the ass. I think it will be *slightly* easier than the status
quo, which requires a giant switch statement depending on the algorithm
to do the same thing.

More concretely: I've been bitten once by a typo in handling/"parsing"
the unwieldy config object (same config object in RStan and PyStan) and
I think using functions would make those kinds of errors much easier to
spot.

Bref, I don't see how it will make things worse and it might make things
better.


Thanks,

Allen

[Bob Carpenter]

I have sympathy --- my problem with the hiearchical structure
is a conceptual one, too.

There are two things I'm suggesting doing differently:

* flattening arguments
* separating top-level commands

My main motivation is to simplify the functions from
the user's perspective by removing argument interactions and
making the functions or command-line invocations more idiomatic.
And also make backward compatibility easier.

So the reason plain-old Euclidean HMC is split from
NUTS is that they take two conflicting config params
(integration time and max depth; whereas the metric is just
a free choice).

The downside is that it increases the number of functions.
This would be bad if a lot of them were used and you had
to constantly choose which one to use.

I can live with stan_mcmc_... and stan_opt_...

I think the full set of top-level commands now would be:

SHORT LONG

stan_lbfgs stan_opt_lbfgs
stan_bfgs stan_opt_bfgs

stan_hmc stan_mcmc_hmc_static
stan_nuts stan_mcmc_hmc_nuts

stan_diagnose stan_diagnose

I think the following are in our short-term future:

stan_bbvb stan_approx_bbvb

stan_rhmc stan_mcmc_hmc_static_riemannian
stan_geonuts stan_mcmc_hmc_nuts_riemannian

and maybe

stan_stochastic_bbvb stan_approx_stochastic_bbvb

stan_mml stan_opt_mml

I'm particularly worried about backward compatibility in
the hierarchical formulation.

I think Andrew/Ben/Rob want to include something to parse
the lme4/glm syntax to implement somehting like:

stan_glm ???
stan_lme ???

- Bob

[Michael Betancourt]

No, no, no, no, no, no. At no point should “stan_hmc”
be used for any command. It’s too general and will
confound the already bad naming conventions in the
literature. Especially when it’s compared to “stan_nuts”
which is not a proper dichotomy.

The more appropriate names would be

stan_mcmc_ehmc_static
stan_mcmc_ehmc_exhaustive

which would then generalize to

stan_mcmc_rhmc_static
stan_mcmc_rhmc_exhaustive

Short names don’t help — if we want a default that
users can call without knowing what they’re doing
then we should just have a “stan_sample” method
that defaults to out current favorite method.

As for glm and lme should these be model format
options? As in you can input a model in the Stan
modeling language, or the glm language, or the
lme language.

[Bob Carpenter]

> On Jun 1, 2015, at 4:06 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> No, no, no, no, no, no. At no point should “stan_hmc”
> be used for any command. It’s too general and will
> confound the already bad naming conventions in the
> literature. Especially when it’s compared to “stan_nuts”
> which is not a proper dichotomy.
>
> The more appropriate names would be
>
> stan_mcmc_ehmc_static

I think the "e" and "static" are defaults and "mcmc" is implied by "hmc",
so it's back to my preferred command,

stan_hmc

by the usual function naming principles (where we treat the noun
"hmc" as an infinitive verb, as in "to verb a noun").

> stan_mcmc_ehmc_exhaustive

Why "exhaustive" for NUTS?

Language is just messy, though we can control these usages if you
really think it's necessary (I don't here and think the longer names
will just be more confusing).

It's completely typical to have a generic term that's also used
for a specific instance. For instance, "blue" might be either a cover
or a contrast to "sky blue" and "navy blue".

And you don't see commands in languages like R that say things like:

mypackagename_linear_regression_nonheteroscedastic_normal_noise(...)

you see commands like lm(), where the lack of "g" means it's plain
old linear regression and the normal noise is taken for granted.

> which would then generalize to
>
> stan_mcmc_rhmc_static
> stan_mcmc_rhmc_exhaustive
>
> Short names don’t help — if we want a default that
> users can call without knowing what they’re doing
> then we should just have a “stan_sample” method
> that defaults to out current favorite method.
>
> As for glm and lme should these be model format
> options? As in you can input a model in the Stan
> modeling language, or the glm language, or the
> lme language.

You'll have to ask Andrew et al. I think the goal was
to have a drop-in replacement for R's glm(), which would
be different than a model-specification option for a
Stan call.

- Bob

[Michael Betancourt]

> I think the "e" and "static" are defaults and "mcmc" is implied by "hmc",
> so it's back to my preferred command,
>
> stan_hmc

I am warning everyone that I will fight this hard. If we want to
abstract to a default then just go with stan_sample or stan_mcmc.
stan_hmc just reinforces this false idea that people have that
HMC is just Euclidean HMC with a static integration time and
every other config is some fancy new algorithm. But it’s not —
HMC is a very general construction and we need to emphasize
this so that people understand we’re just not offering a new
heuristic algorithm but an entirely new perspective of MCMC.

>> stan_mcmc_ehmc_exhaustive
>
> Why "exhaustive" for NUTS?

It’s what were calling a new class of algorithms with dynamic
integration times (as in the integration time expands until
the trajectory “exhausts” exploration on the level set). NUTS
is an example — the binary tree to preserve detailed balance
will be a general element of these methods, although the NUTS
criterion itself may likely change to something more theoretically
motivated.

> And you don't see commands in languages like R that say things like:
>
> mypackagename_linear_regression_nonheteroscedastic_normal_noise(...)
>
> you see commands like lm(), where the lack of "g" means it's plain
> old linear regression and the normal noise is taken for granted.

And people also use “Metropolis” to refer to Random Walk Metropolis
with an isotropic Gaussian proposal, which causes no end of confusion.
For example, I keep hearing “Metropolis is random walk, so all Metropolis
algorithms are bad, right?” which confuses the hell out of people when
you try to respond with HMC being implemented with a Metropolis correction
in practice. Notation and conventions are terrible in statistical software and
that’s no excuse for us to use poor names, too.

Again, if you want to have a default sampling config fine. But don’t make
a default HMC config and call it “stan_hmc” or we’ll be digging out of that
convention for a long, long time.

[Ben Goodrich]

On Monday, June 1, 2015 at 5:24:26 PM UTC-4, Michael Betancourt wrote:

Again, if you want to have a default sampling config fine.  But don’t make
a default HMC config and call it “stan_hmc” or we’ll be digging out of that
convention for a long, long time.

If we make the various configs into formal classes, we can just have the sampling() (or whatever we rename it to) dispatch on the config object in addition to the model object). Then we wouldn't have a proliferation of function names, just a proliferation of class names for the config objects.

Ben

[Ben Goodrich]

I quickly reimplemented the current configuration as S4 classes in R (attached). After a source(), you could do

new("stan_hmc")

or whatever to instantiate an object that has the defaults for HMC. This retains the hierarchical structure of CmdStan calls, validates the input, doesn't have slots that are inapplicable given other options, and does not require us to write too many functions that are callable by the user. Instead, we can write "methods" for the, e.g., sampling() generic function that dispatch according to the class of the config object.

Ben

[Ben Goodrich]

Just going through some of the details of this ...

On Thursday, May 21, 2015 at 6:06:47 PM UTC-4, Bob Carpenter wrote:

1.  Compilation
------------------------------------------------------------
The compilation function compiles a Stan program to a dynamic shared
object (DSO) and optionally saves it.  

compiled_program
rtan::stan_compile(                # DEFAULTS
  source = (<string> || <path>)    # no default
  name = <string>,                 # empty
  boost = <path>,                  # use BH
  eigen = <path>,                  # use RcppEigen
  opt = ("fast" || "debug"),       # "fast"
  save_dso = <bool>,               # FALSE
  debug_code = <bool>              # TRUE
); 

* Is there any reason to keep the name option?

I don't think it hurts anything and the name comes up in some of the printed output.
 

* How long-lived is the "save_dso" option?
    - Where is it saved?

It is kind of confusing. The save_dso option currently just makes the dso part of the stanfit object but you have to manually save the stanfit object to the disk to reuse it in a different session. But this is partially obsoleted by the 2.7.0 behavior if the auto_write option is enabled, in which case it is saving the compiled model object to the disk (in the same directory as the .stan file).
 

* If we allow optimize settings, is there a way to have it scope
  over just the call and not affect the R global build environment?

The inline package doesn't permit passing opt and debug_code to CXX. That has to come from a ~./R/Makevars file.
 

* Is there a way we can have make-like behavior here where if there's
  a saved version and the model file or code have not changed, don't
  do anything?

Already implemented for 2.7.0.
 

* Do we need the other existing compiler options ("presentation1",
  "presentation2", and "small")?

Can't pass them anyway but they were just aliases for -O1, -O2, and -Os.

* Should the save_dso provide a path rather than a boolean?

No. In 2.7.0, it only looks for a foo.rda file in the same directory as foo.stan.
 

* Is it possible to save the DSO across sessions?

Not for the cxxdso object per se but yes for the stanmodel object that has a slot for the cxxdso object.
 

* Should the save_dso defalt be true and if so, use temp?

Yes, that is what it does currently. I don't see much reason to ever specify save_dso = FALSE.
 

* Do we ever need NDEBUG=false?  That is, shouldn't we just always
  have -DNDEBUG set to turn off debugging?

Can't pass it, but in the future I would like NDEBUG to turn off some checks. But that is a Stan thing rather than an interface thing.
 

compiled_program
--------------------
This should only include the compiled Stan program and methods
to do variable conversion, extract names, compute log densities, etc.

Okay, but this would be backwards uncompatible since log_prob() and stuff are currently defined for stanfit objects because they need the data associated with them. So they are going to have to take data arguments. And we should have log_hessian() and maybe laplace_approximation().
 

2. NUTS
------------------------------------------------------------

mcmc_fit
rstan::stan_nuts(

I disagree with having so many suffixed functions for users. I think we should instead code something like

setMethod("sampling", signature(program = "compiled_program", config = "stan_sample"),
          definition = function(program, config, data, pars, init, ...) {
  # implement
})

where an object of stan_nuts (or stan_fixed or whatever) class could be passed to the config argument or if not, a default one would be instantiated and used.
 

* Return value just includes the result in such a way that it
  can be extracted as before.

OK
 

* Return value won't include the compiled program (or source).

I don't see how including the compiled_program would hurt anything
 

* Should return value include timing information?

Some people have asked for it, although currdently we don't do so programatically
 

* There's no append_samples = <bool>

OK
 

* There's no option to compile the model from within the stan_nuts
  call.

I don't think there should be a stan_nuts() call, but either way, it should pass the compiled_model in.
 

* data defaults to top-level environment with lazy read (not how
  it's done now);  list is list of values as before

I'll implement that for 2.7.0 since no one objected to what I posted before. It is technically the environments in R's lexical order rather than always the top-level environment.
 

* Instead of "chains" and a start ID, there's just a "chain_ids"

Currently (and I think it should stay this way in the future) sampling() only does one chain, so it only needs to know the singular chain_id for the chain it is doing. Whatever is driving it whether it is in serial or parallel should be responsible for passing it the seed, chain_id, etc.
 

* Default num_warmup and num_samples very small by default (but shouldn't
  be so small as to trigger scary messages about adaptation)

Are you proposing changing the defaults of num_warmup and num_samples (1000 each)? If so, then we should do that for all the interfaces.
 

* Metric is specified by name of type or by value;  the names default
  to unit values for init
* Those adapt parameters are all for adapt_stepsize, right?

Yeah, although I think the various configs should have their own class structure.
 

* Are there adaptation parameters for mass_matrix, like regularization?

Not currently but we could add them.
 

* params_to_save defaults to NA (or all), or a positive list of whole
  parameter names.

I'm not sure whether it is worth renaming the pars argument params_to_save, in part because in 2.7.0 there is an include = TRUE argument, so you could specify pars = "Sigma", include = FALSE to exclude Sigma but include everything else.
 

* sample_file gets a copy of the output written to it, but return
  object still created --- should this just move so that the return
  object has a save on it?  Or should setting this cause the return of
  this function to be NA?

I think it should still return the regular stuff, in addition to writing the sample_file(s).
 

mcmc_fit
--------------------


3. EHMC
------------------------------------------------------------
Plain old HMC looks just like NUTS, but with max_treedepth replaced
by integration time:

I would say they should be siblings in the class hierarchy for stan_config objects but get called the same way.
 

4. Optimize
------------------------------------------------------------

optimization_fit
stan_lbfgs(

Same deal. I think stan_lbfgs should be a class that inherits from stan_bfgs class with an additional slot for history_size but the user would call something like

stan_optimize(compiled_program, config, data, init, verbose, ...)

* Should verbose be a boolean?

Yes
 

* What does current as_vector do? Current doc says "flag indicating whether a
  vector is used to for the point estimate found. A list can be
  used instead by specifying it to be FALSE"

If as_vector = TRUE, then only a vector of optimal parameters is returned. If as_vector = FALSE, then there is more stuff returned with the list.
 

* What's the behavior of tolerances if more than one's specified?  I'd
  think it should be read as AND or as OR.

Don't know.
 

* No option to specify parameter return as vector --- use a compiled
  model's conversion.

But the compiled_model doesn't know what the optimal parameter values are. At least it doesn't currently.
 

optimization_fit
--------------------
    $par : parameters (as list)
    $value : lp__

And $hessian optionally
 

5. Diagnose
------------------------------------------------------------

diagnostics
rstan::stan_diagnose(

I think stan_diagnose should be a subclass of stan_config with epsilon and error_threshold that can be passed to stan_sample() or stan_optimize() like any other stan_config object.

diagnostics
--------------------
?? not sure what should go in here --- it needs to be printable

It doesn't exactly print (in the R sense currently) but the gradients get written to R's stream from C++. I don't think the gradients get saved anywhere currently, although we could do that.

Ben

[Bob Carpenter]

Thanks --- that was really helpful. More inline below.

> On Jun 2, 2015, at 12:47 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> Just going through some of the details of this ...
>
> On Thursday, May 21, 2015 at 6:06:47 PM UTC-4, Bob Carpenter wrote:
> 1. Compilation
> ------------------------------------------------------------
> The compilation function compiles a Stan program to a dynamic shared
> object (DSO) and optionally saves it.
>
> compiled_program
> rtan::stan_compile( # DEFAULTS
> source = (<string> || <path>) # no default
> name = <string>, # empty
> boost = <path>, # use BH
> eigen = <path>, # use RcppEigen
> opt = ("fast" || "debug"), # "fast"
> save_dso = <bool>, # FALSE
> debug_code = <bool> # TRUE
> );
> * Is there any reason to keep the name option?
>
> I don't think it hurts anything and the name comes up in some of the printed output.

What it hurts is readability of the doc and maintainability of
the code.

> * How long-lived is the "save_dso" option?
> - Where is it saved?
>
> It is kind of confusing. The save_dso option currently just makes the dso part of the stanfit object but you have to manually save the stanfit object to the disk to reuse it in a different session.

I'd like to make the more fine-grained objects available, like
being able to take the DSO and give it new data without recompiling.
(More on this below.)

> But this is partially obsoleted by the 2.7.0 behavior if the auto_write option is enabled, in which case it is saving the compiled model object to the disk (in the same directory as the .stan file).
>
> * If we allow optimize settings, is there a way to have it scope
> over just the call and not affect the R global build environment?
>
> The inline package doesn't permit passing opt and debug_code to CXX. That has to come from a ~./R/Makevars file.
>
> * Is there a way we can have make-like behavior here where if there's
> a saved version and the model file or code have not changed, don't
> do anything?
>
> Already implemented for 2.7.0.

Fantastic. It's working now? I thought you'd said it was very fragile.

> * Do we need the other existing compiler options ("presentation1",
> "presentation2", and "small")?
>
> Can't pass them anyway but they were just aliases for -O1, -O2, and -Os.

What does "can't pass them" mean? If they're impossible, we should remove
them from the doc.

> * Should the save_dso provide a path rather than a boolean?
>
> No. In 2.7.0, it only looks for a foo.rda file in the same directory as foo.stan.

OK

> * Is it possible to save the DSO across sessions?
>
> Not for the cxxdso object per se but yes for the stanmodel object that has a slot for the cxxdso object.

Would it be possible? It would be super great if we had
a DSO wrapper object and it was saveable, and an instantiated
object that could be saved.

> * Should the save_dso defalt be true and if so, use temp?
>
> Yes, that is what it does currently. I don't see much reason to ever specify save_dso = FALSE.

OK, then I think we should get rid of the option.

> * Do we ever need NDEBUG=false? That is, shouldn't we just always
> have -DNDEBUG set to turn off debugging?
>
> Can't pass it, but in the future I would like NDEBUG to turn off some checks. But that is a Stan thing rather than an interface thing.

Couldn't we do a preliminary call to set_cppo() or the like, or
would that run afoul of CRAN policy on explicitness of side effects?

> compiled_program
> --------------------
> This should only include the compiled Stan program and methods
> to do variable conversion, extract names, compute log densities, etc.
>
> Okay, but this would be backwards uncompatible since log_prob() and stuff are currently defined for stanfit objects because they need the data associated with them. So they are going to have to take data arguments. And we should have log_hessian() and maybe laplace_approximation().

I would very much like to move away from packing so much into
the fit object. I think we should maintain backward compatibility
until Stan 3, but there's no reason we can't introduce things the way
we want them and mark others as deprecated.

You're right that I missed a level here.

1. Compiled program --- this is the DSO thing.

2. Instantiated program class --- this is the program wedded
to its data.

We actually need the instance in (2) to do anything I suggested.

> 2. NUTS
> ------------------------------------------------------------
>
> mcmc_fit
> rstan::stan_nuts(
>
> I disagree with having so many suffixed functions for users. I think we should instead code something like
>
> setMethod("sampling", signature(program = "compiled_program", config = "stan_sample"),
> definition = function(program, config, data, pars, init, ...) {
> # implement
> })
>
> where an object of stan_nuts (or stan_fixed or whatever) class could be passed to the config argument or if not, a default one would be instantiated and used.

OK --- Rob unwound what all this meant for me, which I think means
that it'll be used as follows:

> program <- stan_model(...);
> fit <- sampling(program, config, ...);

where the config is a required argument to sampling().

Then the idea's to do what Michael likes and have a big config object
that works for everything --- optimization, MCMC, diagnosis, and whatever
else comes up? Or do you like breaking things out into sampling(),
optimizing(), diagnosing(), vbing() etc.?

I'd want to make things flatter, but I can live with this.

With this, would we need the top-level stan() object to make
people happy (again, I don't want to rip it out ASAP --- not until
an RStan 3 was solidified).

> * Return value just includes the result in such a way that it
> can be extracted as before.
>
> OK
>
> * Return value won't include the compiled program (or source).
>
> I don't see how including the compiled_program would hurt anything

It makes the objects more complicated and it makes it hard to
build a valid return object without one. Then functions need to all
be defined conditionally onw hether one's there or not if it's optional,
which is hard to debug and maintain and doc.

> * Should return value include timing information?
>
> Some people have asked for it, although currdently we don't do so programatically

This isn't a big deal.

> * There's no append_samples = <bool>
>
> OK
>
> * There's no option to compile the model from within the stan_nuts
> call.
>
> I don't think there should be a stan_nuts() call, but either way, it should pass the compiled_model in.

OK.

> * data defaults to top-level environment with lazy read (not how
> it's done now); list is list of values as before
>
> I'll implement that for 2.7.0 since no one objected to what I posted before. It is technically the environments in R's lexical order rather than always the top-level environment.

Sorry -- I didn't mean top-level environment, I meant whatever
environment the RStan method/function was called from.

> * Instead of "chains" and a start ID, there's just a "chain_ids"
>
> Currently (and I think it should stay this way in the future) sampling() only does one chain, so it only needs to know the singular chain_id for the chain it is doing. Whatever is driving it whether it is in serial or parallel should be responsible for passing it the seed, chain_id, etc.

I'm guessing users like Andrew (or even me) are going to want something
that runs multiple chains and returns a single object on which we can print(),
plot(), and extract(). So I think this has to deal with multiple chains,
or at least we need top-level drivers that are included that do that.

> * Default num_warmup and num_samples very small by default (but shouldn't
> be so small as to trigger scary messages about adaptation)
>
> Are you proposing changing the defaults of num_warmup and num_samples (1000 each)? If so, then we should do that for all the interfaces.

Yes, I think we agreed on 100 for all the interfaces.
We should make sure everyone agrees. Too many people turn on 2000 and
then go away for the weekened.

> * Metric is specified by name of type or by value; the names default
> to unit values for init
> * Those adapt parameters are all for adapt_stepsize, right?
>
> Yeah, although I think the various configs should have their own class structure.

I agree.

> * Are there adaptation parameters for mass_matrix, like regularization?
>
> Not currently but we could add them.

I think we'd need to add them through the Stan C++ side
first.

Some of these questions were intended more generally, not about
the current state of RStan or 2.7.

> * params_to_save defaults to NA (or all), or a positive list of whole
> parameter names.
>
> I'm not sure whether it is worth renaming the pars argument params_to_save, in part because in 2.7.0 there is an include = TRUE argument, so you could specify pars = "Sigma", include = FALSE to exclude Sigma but include everything else.

OK.

> * sample_file gets a copy of the output written to it, but return
> object still created --- should this just move so that the return
> object has a save on it? Or should setting this cause the return of
> this function to be NA?
>
> I think it should still return the regular stuff, in addition to writing the sample_file(s).
>
> mcmc_fit
> --------------------
>
>
> 3. EHMC
> ------------------------------------------------------------
> Plain old HMC looks just like NUTS, but with max_treedepth replaced
> by integration time:
>
> I would say they should be siblings in the class hierarchy for stan_config objects but get called the same way.

Looks like Michael strongly agrees with this. We absolutely
have to get Michael's approval on whatever we call anything that
smacks of "HMC". I'm very confused as to why Michael thinks
everyone's confused, but that's the nature of confusion.

> 4. Optimize
> ------------------------------------------------------------
>
> optimization_fit
> stan_lbfgs(
>
> Same deal. I think stan_lbfgs should be a class that inherits from stan_bfgs class with an additional slot for history_size but the user would call something like
>
> stan_optimize(compiled_program, config, data, init, verbose, ...)

OK --- I don't know what stan_bfgs is. But in general, I don't like
the idea of lbfgs inheriting from bfgs because it's not an ISA relationship.
They could both inherit from a shared class that does what's in common, which
is how Michael has a lot of the MCMC classes set up in Stan C++.

> * Should verbose be a boolean?
>
> Yes
>
> * What does current as_vector do? Current doc says "flag indicating whether a
> vector is used to for the point estimate found. A list can be
> used instead by specifying it to be FALSE"
>
> If as_vector = TRUE, then only a vector of optimal parameters is returned. If as_vector = FALSE, then there is more stuff returned with the list.

Do we ever want as_vector = TRUE? If not, let's get rid of it.
That is, if it's not too heavy, I think it's easier to return the whole
thing the same way all the time, then have an extractor or converter
for the result.

> * What's the behavior of tolerances if more than one's specified? I'd
> think it should be read as AND or as OR.
>
> Don't know.

This just needs to go in the doc. It's explained in either the Stan
manual or the CmdStan doc, but I can never remember.

> * No option to specify parameter return as vector --- use a compiled
> model's conversion.
>
> But the compiled_model doesn't know what the optimal parameter values are. At least it doesn't currently.

Right, those would have to be passed in.

> optimization_fit
> --------------------
> $par : parameters (as list)
> $value : lp__
>
> And $hessian optionally
>
> 5. Diagnose
> ------------------------------------------------------------
>
> diagnostics
> rstan::stan_diagnose(
>
> I think stan_diagnose should be a subclass of stan_config with epsilon and error_threshold that can be passed to stan_sample() or stan_optimize() like any other stan_config object.

Unless that subclass is a common subclass, then I don't like it
for the reasons above. Only subclass for ISA relationships.

> diagnostics
> --------------------
> ?? not sure what should go in here --- it needs to be printable
>
> It doesn't exactly print (in the R sense currently) but the gradients get written to R's stream from C++. I don't think the gradients get saved anywhere currently, although we could do that.

If it returns an object, that object should have a print() method
and ideally a plot() method. Then if you just call it without assigning,
it'll do the right thing and print, right?

- Bob

[Allen B. Riddell]

It seems like we agree on most things.

Discussion seems to be focused on:

- fragmentation/coagulation of configuration object and/or functions
- how thin should the stan_fit object be
- does "name" survive
- behavior of sample_file

We might keep in mind that we can probably accommodate some differences
between interfaces. For example, on `sample_file` I don't know of a
Python package that both returns values and writes out the same values
to a csv file in the same step (esp. in a world that's moving away from
CSV files and towards things like JSON and HDF5), but it probably
doesn't matter if one interface does something slightly different here.

(a few comments inline below)

Perhaps we can consolidate these on the Wiki? Including the points to be
resolved?

Best,

Allen

On 06/02, Ben Goodrich wrote:
> * Is there any reason to keep the name option?
> >
>
> I don't think it hurts anything and the name comes up in some of the
> printed output.
>
>

I'd vote emphatically for removing `name`. Models aren't generally named
in R or Python -- we don't write lm0 <- lm(y ~ x, name="mymodel"). If
the user wants to name their model, they are free to do so. Insofar as
Stan is concerned, the identity of the model is the compiled model's C++
code (or the hash thereof).

> > * How long-lived is the "save_dso" option?
> > - Where is it saved?
> >
>
> It is kind of confusing. The save_dso option currently just makes the dso
> part of the stanfit object but you have to manually save the stanfit object
> to the disk to reuse it in a different session. But this is partially
> obsoleted by the 2.7.0 behavior if the auto_write option is enabled, in
> which case it is saving the compiled model object to the disk (in the same
> directory as the .stan file).
>

Is auto_write (aside: is this the best variable name?) going to be "on"
by default (yes!).

> ... <snip> ...

>
> > * sample_file gets a copy of the output written to it, but return
> > object still created --- should this just move so that the return
> > object has a save on it? Or should setting this cause the return of
> > this function to be NA?
> >
>
> I think it should still return the regular stuff, in addition to writing
> the sample_file(s).
>

Is there precedent for this in any R packages? I can't think of any
Python packages that do both.

[Ben Goodrich]

Okay, a few meta comments:

• The set_cppo() function (and doc) are currently gone because of the considerable probability of interfering with other R packages a user might install. I guess I should bring set_cppo() (and doc) back but only have it throw a warning that it doesn't do anything any more. We can put on the wiki instructions on how to modify ~/.R/Makevars but there is not much of a difference between -O3 and the default of -O2.
• We can presumably create a class that combines the compiled code with R objects that correspond to the declarations in the data block of the .stan program. But currently I don't know exactly how to implement something like that in R code. My guess is that the code we have now does it, but in a transitory way since it is not a stand-alone thing.
• I put a rough implementation of my "config S4 objects" idea in the other post https://groups.google.com/d/msg/stan-dev/jX-D619WW0o/K_okbvXaBKgJ . Those objects are not any bigger than they have to be in the sense that e.g. stan_bfgs only has BFGS-related slots but not num_warmup, etc. that are only relevant for MCMC.
• If we internally need to pad a config object so that it has every possible configuration detail (even if the irrelevant ones are not accessed on the C++ side), then we could do that without involving the user.
• About the ISA thing, it depends on what the definition of is is, and R defines it a bit differently than Java does. In this sense, a stan_lbfgs class "is" a stan_bfgs class because all of the slots of stan_bfgs are applicable to stan_lbfgs but stan_lbfgs adds one additional slot for history_size. So, the logic is more like "the child class includes everything that the parent class has, plus more" than a "the child class is a special case of the parent class". This does allow for minimally code validation because a stan_lbfgs has to be a valid stan_bfgs in addition to having a positive history_size. I think you are saying there should be a stan_bfgs_or_lbfgs class that is a parent of each, but then an object of stan_bfgs_or_lbfgs class can be created even though it is insufficient to know what to do with it (i.e. call the BFGS algorithm or the LBFGS algorithm).
  

Details ...

On Tuesday, June 2, 2015 at 2:01:41 PM UTC-4, Bob Carpenter wrote:

> * Is there any reason to keep the name option?
>
> I don't think it hurts anything and the name comes up in some of the printed output.

What it hurts is readability of the doc and maintainability of
the code.

I think those are minor relative to being to the print refering to "my_cool_model" rather than a md5 hash.
 

> * Is it possible to save the DSO across sessions?
>
> Not for the cxxdso object per se but yes for the stanmodel object that has a slot for the cxxdso object.

Would it be possible?  It would be super great if we had
a DSO wrapper object and it was saveable, and an instantiated
object that could be saved.

Maybe. I don't see much utility of saving the instantiated object with the data to the disk, since your data are presumably already saved to the disk along with the script that loads them. The good thing about the foo.rda files is that they are of negligible size on the disk.
 

Then the idea's to do what Michael likes and have a big config object
that works for everything --- optimization, MCMC, diagnosis, and whatever
else comes up?  Or do you like breaking things out into sampling(),
optimizing(), diagnosing(), vbing() etc.?  

I'd want to make things flatter, but I can live with this.

I think there should be sampling() and optimizing() --- perhaps with reverbing --- that would take each take a subset of stan_config classes. To my mind, VB is optimizing and just needs a (future) stan_vb config to distinguish it from ML. Both sampling() and optimizing() could accept a stan_diagnose config since at the moment it just streams the gradient comparison, but we could do something different if we ever had algorithm-specific diagnostics. Then stan_hmc, stan_rmhmc, etc. would be valid configs for sampling(). And if Michael needs us to utilize an uber-class internally to call the C++, then fine (this detail has never struck me as a big concern).
 

With this, would we need the top-level stan() object to make
people happy (again, I don't want to rip it out ASAP --- not until
an RStan 3 was solidified).

Yeah, stan() would basically look the same but calling sampling() with config=new("stan_hmc") as many times as there are chains (optionally in parallel). I think it is no more fragile than Andrew's stan() function on mc-stan.org.
 

> * Default num_warmup and num_samples very small by default (but shouldn't
>   be so small as to trigger scary messages about adaptation)
>
> Are you proposing changing the defaults of num_warmup and num_samples (1000 each)? If so, then we should do that for all the interfaces.

Yes, I think we agreed on 100 for all the interfaces.
We should make sure everyone agrees.  Too many people turn on 2000 and
then go away for the weekened.

I agree that it should be the same for all interfaces. I would never call it with 100 because that is too few to do the convergence diagnostics I prefer. But more broadly, I disagree with the goal of making default arguments be optimal; default arguments should be on the safe side. If someone knows that they are doing MCMC in a situation where ML would work just fine, then they can tweak the options to reduce the runtime.
 

> Plain old HMC looks just like NUTS, but with max_treedepth replaced
> by integration time:
>
> I would say they should be siblings in the class hierarchy for stan_config objects but get called the same way.

Looks like Michael strongly agrees with this.  We absolutely
have to get Michael's approval on whatever we call anything that
smacks of "HMC".  I'm very confused as to why Michael thinks
everyone's confused, but that's the nature of confusion.

In my initial implementation, I made them be almost exactly like CmdStan, but if we are renaming CmdStan arguments to then we can rename them on the interface side.

> * What does current as_vector do? Current doc says "flag indicating whether a
>   vector is used to for the point estimate found. A list can be
>   used instead by specifying it to be FALSE"
>
> If as_vector = TRUE, then only a vector of optimal parameters is returned. If as_vector = FALSE, then there is more stuff returned with the list.

Do we ever want as_vector = TRUE?  If not, let's get rid of it.
That is, if it's not too heavy, I think it's easier to return the whole
thing the same way all the time, then have an extractor or converter
for the result.

If I were writing it from scratch, I wouldn't include an as_vector flag, but I'm not sure it is worth breaking any R code that might be using it.

> diagnostics
> --------------------
> ?? not sure what should go in here --- it needs to be printable
>
> It doesn't exactly print (in the R sense currently) but the gradients get written to R's stream from C++. I don't think the gradients get saved anywhere currently, although we could do that.

If it returns an object, that object should have a print() method
and ideally a plot() method.  Then if you just call it without assigning,
it'll do the right thing and print, right?

Currently, I don't think C++ returns anything useful. It might be better to not stream but return the gradient comparison as a matrix that the user would have to print().

Ben

[Ben Goodrich]

On Tuesday, June 2, 2015 at 3:19:14 PM UTC-4, Allen B. Riddell wrote:

Is auto_write (aside: is this the best variable name?) going to be "on"
by default (yes!).

No. CRAN would consider it malicious if functions were silently writing to the disk by default.
 

> > * sample_file gets a copy of the output written to it, but return
> >   object still created --- should this just move so that the return
> >   object has a save on it?  Or should setting this cause the return of
> >   this function to be NA?
> >
>
> I think it should still return the regular stuff, in addition to writing
> the sample_file(s).
>

Is there precedent for this in any R packages? I can't think of any
Python packages that do both.

None that come to mind. The problem is that read_stan_csv() makes a more limited object than does stan() currently. I usually use the sample_file option when doing more time consuming runs so that I can look at the draws in another process but (assuming it is going well) I still want the fully functional stanfit object whenever it is done.

Ben

[Allen B. Riddell]

On 06/02, Ben Goodrich wrote:

> On Tuesday, June 2, 2015 at 3:19:14 PM UTC-4, Allen B. Riddell wrote:
> >
> > Is auto_write (aside: is this the best variable name?) going to be "on"
> > by default (yes!).
> >
>
> No. CRAN would consider it malicious if functions were silently writing to
> the disk by default.
>

Couldn't one use the user-specific temporary files directory. I think
Windows, OS X, and Linux all have one. On Linux it's ~/.cache and on
Windows it is C:\users\username\AppData\Local\Temp

I'm strongly in favor of creating a LRU cache of compiled models in that
directory by default while making sure not to use more than N megabytes,
where N is the same as used by Chrome/Firefox by default for its
temporary files cache.

I think this would be a great feature and help us with some anecdotal
speed comparisons vis-a-vis other programs. I feel like people (myself
included) often are recompiling identical models because they forgot to
save them.

Best,

Allen

[Allen B. Riddell]

On 06/02, Ben Goodrich wrote:

So it seems like there's another issue here we should discuss -- i.e.,
standard/recommended way to serializing and deserializing stan_model and
stan_fit objects. Since if you need to get a full stan_fit object (or
stan_model object) there should definitely be a way to do that!

On the more narrow point about potentially getting rid of sample_file --
I think the proposed change would have the user in R or Python take
responsibility for writing the data in the stan_fit object to a csv file
after sampling/optimizing. We could provide a helper function to assist
the user -- but the rationale for the change would be that this writing
should be a separate step and not something done as part of the fitting.

Best,

Allen

[Bob Carpenter]

> On Jun 2, 2015, at 3:44 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> Okay, a few meta comments:
> • The set_cppo() function (and doc) are currently gone because of the considerable probability of interfering with other R packages a user might install. I guess I should bring set_cppo() (and doc) back but only have it throw a warning that it doesn't do anything any more. We can put on the wiki instructions on how to modify ~/.R/Makevars but there is not much of a difference between -O3 and the default of -O2.

OK. We should test the performance of -O2 vs. -O3.

But I agree that it's good not to stomp on other people's config
through a function if it's not clear to the user what they're doing.

> • We can presumably create a class that combines the compiled code with R objects that correspond to the declarations in the data block of the .stan program. But currently I don't know exactly how to implement something like that in R code. My guess is that the code we have now does it, but in a transitory way since it is not a stand-alone thing.

I don't get the context for this. I just want to be able
to grab:

1. compiled class for Stan program

2. class for Stan program instantiated with data

We must be doing at least (2) now because the functionality's there.
Or maybe it just gets reconsutructed each time?

> • I put a rough implementation of my "config S4 objects" idea in the other post https://groups.google.com/d/msg/stan-dev/jX-D619WW0o/K_okbvXaBKgJ . Those objects are not any bigger than they have to be in the sense that e.g. stan_bfgs only has BFGS-related slots but not num_warmup, etc. that are only relevant for MCMC.

Understood.

> • If we internally need to pad a config object so that it has every possible configuration detail (even if the irrelevant ones are not accessed on the C++ side), then we could do that without involving the user.

> • About the ISA thing, it depends on what the definition of is is,

The relevant definition of "is" is that "A is B" if and only if whatever
you can do with a B you can do with an A. So it's not only having the same
slots, but having them mean the same thing. That is, if you have A extend B,
then you should be able to run A as a B object to do everything B does and you
should be able to assign an A to a B reference.

> and R defines it a bit differently than Java does.

Java doesn't legislate this in its design and lots of people abuse
the inheritance hierarchy to do things that should be done with
encapulation. That is, you see lots of people extend the library String
class to do things like add a number, so you get a String with a number.
That's not technically wrong if you think of your extension class as a String
with a number, but that's not usually what you want.

> In this sense, a stan_lbfgs class "is" a stan_bfgs class because all of the slots of stan_bfgs are applicable to stan_lbfgs but stan_lbfgs adds one additional slot for history_size. So, the logic is more like "the child class includes everything that the parent class has, plus more" than a "the child class is a special case of the parent class". This does allow for minimally code validation because a stan_lbfgs has to be a valid stan_bfgs in addition to having a positive history_size. I think you are saying there should be a stan_bfgs_or_lbfgs class that is a parent of each, but then an object of stan_bfgs_or_lbfgs class can be created even though it is insufficient to know what to do with it (i.e. call the BFGS algorithm or the LBFGS algorithm).

So there's no way to make a class purely abstract (and hence uninstantiatable)
as in Java or C++?

- Bob

[Michael Betancourt]

Dammit, we need to have as detailed conversations while
I’m conferencing. I’ll try to look through the thread in
depth tomorrow during the talks, but I’ll note that I’ve
always liked classes with default config. So if we
had a Stan class it would get initialized with a valid
config and then we’d have mutators that checked
validity ensuring that the internal config is always
valid. That said, this does go against the
static-config-at-construction pattern we mostly use
internally to Stan and I’m not sure how we’d design
a parallel pattern for CmdStan.

[Bob Carpenter]

> On Jun 2, 2015, at 5:14 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> Dammit, we need to have as detailed conversations while
> I’m conferencing.

I think it's hard to have these detailed discussions on our
low-bandwidth conference calls with so many people online.
I often just lose whole bits of the conversation when I'm online
because people are talking over whoever's speaking in the
background, which is fine in Andrew's office, but garbles everything
online.

> I’ll try to look through the thread in
> depth tomorrow during the talks, but I’ll note that I’ve
> always liked classes with default config. So if we
> had a Stan class it would get initialized with a valid
> config and then we’d have mutators that checked
> validity ensuring that the internal config is always
> valid. That said, this does go against the
> static-config-at-construction pattern we mostly use
> internally to Stan and I’m not sure how we’d design
> a parallel pattern for CmdStan.

The reason to make classes immutable is that it's easier to
guarantee they're always in safe and usable states. And it
automatically makes them threadsafe, which would otherwise require
synchronization on the gets/sets.

It's often hard to do with getters/setters if you have linked
states to get and set. For instance, if you have a class

struct {
real l_;
real u_;
}

with a constraint that l_ < u_, then you can't update that without
setting l_ and u_ at the same time because they're not independent
in validation. The more interactions, the harder the getter/setter
design (and the harder the validation in general).

- Bob

[Michael Betancourt]

Right, but a hierarchical config wouldn’t have any interactions

per say.  You could modify NUTS parameters even if NUTS

wasn’t the active HMC implementation.  In other words, you

have a big ol’ tree of parameters and the mutators would

ensure that each node has a valid value (validity being

defined by each parameter definition) and then at runtime

the config gets executed and selects out a single subtree

of active parameters.  Because the global tree is valid

each subtree will be valid.

Anyways, I’m just saying that a global config object

would be easy to manipulate safely.  It doesn’t get

around the overhead for learning a hierarchical

structure.  Also, there’s the interface consistency.

I think if we went with something like this we’d

have to not use methods but rather have functions

that take in a config class,

stan_config

stan_config.hmc.stepsize.set(0.0001)

sample(stan_config)

which would then roughly parallel the CmdStan

execution of 

./model sample stan_config_string

[Ben Goodrich]

On Tuesday, June 2, 2015 at 7:41:26 PM UTC-4, Michael Betancourt wrote:

I think if we went with something like this we’d

have to not use methods but rather have functions

that take in a config class,

stan_config

stan_config.hmc.stepsize.set(0.0001)

sample(stan_config)

which would then roughly parallel the CmdStan

execution of 

./model sample stan_config_string

I don't see the case for making the user deal with an object that has every possible configuration option. If the user wants to do HMC, then

config <- new("stan_hmc", stepsize = 0.0001) # is validated in R
print(config) # outputs something helpful but tailored to HMC
posterior <- sampling(model, config)

and then sampling() would internally coerce that config into whatever is needed by C++, including possibly filling in all the extra unused configuration options.

Ben

[Michael Betancourt]

* Should return value include timing information?

Some people have asked for it, although currdently we don't do so programatically

This changes in the writer refactor, where timing information is passed

to the interfaces via writer callbacks.  It will then be trivial for RStan to

store the timings or do whatever it likes.

[Michael Betancourt]

Are you proposing changing the defaults of num_warmup and num_samples (1000 each)? If so, then we should do that for all the interfaces.

Yes, I think we agreed on 100 for all the interfaces.
We should make sure everyone agrees.  Too many people turn on 2000 and

then go away for the weekend.

I think 100 might be too small.  It might be fine for social science models

but for complex likelihoods (especially with ODEs) the adaptation typically

needs longer than 100.  I’d suggest 250 minimum, with a warning to users

that they might need more if they’re seeing poor sampling.

Looks like Michael strongly agrees with this.  We absolutely
have to get Michael's approval on whatever we call anything that
smacks of "HMC".  I'm very confused as to why Michael thinks
everyone's confused, but that's the nature of confusion.

Once again, the problem is that people currently equate “Hamiltonian

Monte Carlo” with what should more formally be called “Metropolis-corrected 

Euclidean Hamiltonian Monte Carlo with static integration time”.  So when

you tweak any of those properties people think it’s a different algorithm

and we lose the concept of Hamiltonian Monte Carlo being this very

general approach.  This is a ubiquitous problem in MCMC  as it also

happens with Metropolis-Hastings (people confusing the general 

Metropolis-Hastings construction vs particular proposals).

I want the naming to emphasize that any default is just _an_ implementation

of Hamiltonian Monte Carlo not _the_ implementation of Hamiltonian Monte

Carlo.  Having a “stan_hmc” option does the latter.

[Michael Betancourt]

We might keep in mind that we can probably accommodate some differences
between interfaces. For example, on `sample_file` I don't know of a
Python package that both returns values and writes out the same values
to a csv file in the same step (esp. in a world that's moving away from
CSV files and towards things like JSON and HDF5), but it probably
doesn't matter if one interface does something slightly different here.

This is another thing that changes in the writer refactor as all

sample output gets passed to the interfaces piece by piece

by the callbacks.

But this does raise an important issue — in the new refactored

architecture each interface will be responsible for creating

writers and var_contexts which may or may not require their

own arguments.  So we can’t talk about a single config object

in isolation, rather we’ll have “internal” arguments and “interface”

arguments and we should keep this dichotomy in mind.

[Michael Betancourt]

I think there should be sampling() and optimizing() --- perhaps with reverbing --- that would take each take a subset of stan_config classes. To my mind, VB is optimizing and just needs a (future) stan_vb config to distinguish it from ML. Both sampling() and optimizing() could accept a stan_diagnose config since at the moment it just streams the gradient comparison, but we could do something different if we ever had algorithm-specific diagnostics. Then stan_hmc, stan_rmhmc, etc. would be valid configs for sampling(). And if Michael needs us to utilize an uber-class internally to call the C++, then fine (this detail has never struck me as a big concern).
 

I disagree — I think ever new algorithmic approach (optimization, MCMC, VB, EP, MML, diagnostics) should

be a separate call.  VB is very much fundamentally different from just an optimization of some objective function.

[Ben Goodrich]

But given the existence of a configuration object that uniquely defines the approach, what is the benefit to users of making them learn to call a different function whose name corresponds to part of the configuration? VB is fundamentally different from penalized ML but how is optimization with VB configuration fundamentally different from optimization with penalized ML configuration?

Ben

[Michael Betancourt]

All of the BBVB tuning parameters, for one.  But eventually also possibly

the variation distribution, or the exact variational algorithm or how the

variation expectations are reported on the constrained scale.

[Ben Goodrich]

On Wednesday, June 3, 2015 at 9:48:02 AM UTC-4, Michael Betancourt wrote:

All of the BBVB tuning parameters, for one.  But eventually also possibly

the variation distribution, or the exact variational algorithm or how the

variation expectations are reported on the constrained scale.

But those are all aspects of the VB config. It would specify all the tuning parameters that are VB-specific. And we could add flags to it as necessary to select the variation distribution or how it is minimized. I think how the expectations are reported would be a print option, but in principle that could be part of the VB config too.

I just don't see the utility of having the user call

config <- new("stan_vb_config")
stan_vb(program, config)
stan_hmc(program, config) # error

when that config already implies the VB approach and is only valid input for a stan_vb() function.

Ben
 

On Jun 3, 2015, at 3:46 PM, Ben Goodrich <goodri...@gmail.com> wrote:

On Wednesday, June 3, 2015 at 2:30:54 AM UTC-4, Michael Betancourt wrote:

I think there should be sampling() and optimizing() --- perhaps with reverbing --- that would take each take a subset of stan_config classes. To my mind, VB is optimizing and just needs a (future) stan_vb config to distinguish it from ML. Both sampling() and optimizing() could accept a stan_diagnose config since at the moment it just streams the gradient comparison, but we could do something different if we ever had algorithm-specific diagnostics. Then stan_hmc, stan_rmhmc, etc. would be valid configs for sampling(). And if Michael needs us to utilize an uber-class internally to call the C++, then fine (this detail has never struck me as a big concern).
 

I disagree — I think ever new algorithmic approach (optimization, MCMC, VB, EP, MML, diagnostics) should

be a separate call.  VB is very much fundamentally different from just an optimization of some objective function.

But given the existence of a configuration object that uniquely defines the approach, what is the benefit to users of making them learn to call a different function whose name corresponds to part of the configuration? VB is fundamentally different from penalized ML but how is optimization with VB configuration fundamentally different from optimization with penalized ML configuration?

Ben

--
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.