ragged array assignment (not declaration)

[bob-carpenter]

From Ben Goodrich on stan-dev:

OK. I invented a new rule (locally): You are allowed to assign anything to a matrix with 0 rows and 0 columns. So, this now works:

model {
  matrix[0,0] X[2];
  vector[4] ones;
  ones[1] <- 1;
  ones[2] <- 1;
  ones[3] <- 1;
  ones[4] <- 1;
  X[1] <- diag_matrix(ones);
  theta ~ normal(0,1);
}

But it would probably be better to support (outside of the data and parameters blocks)

matrix X[2]; 
// parse to vector<Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic> > X(2);
// i.e. no fill value after the 2

rather than making users declare

matrix[0,0] X[2];

I think we're good in that case in the sense that this runs so I guess the column vector gets promoted:

#include <vector>
#include <Eigen/Dense>

int main() {
  std::vector<Eigen::Matrixd> > x(2);
  Eigen::MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  x[0] = m;
  x[1] = m.col(1);
  return 0;
}

However, the comment in src/stan/math/prim/mat/fun/assign.hpp does not seem to match what the code does. I had to do trial and error to figure out that this was being called rather than the one above it that allegedly fires when the dimensions do not match.

/** 
 * Copy the right-hand side's value to the left-hand side 
 * variable. 
 * 
 * The <code>assign()</code> function is overloaded.  This 
 * instance will be called for arguments that are both 
 * <code>Eigen::Matrix</code> types and whose shapes match.  The 
 * shapes are specified in the row and column template parameters. 
 * 
 * @tparam LHS Type of left-hand side matrix elements. 
 * @tparam RHS Type of right-hand side matrix elements. 
 * @tparam R Row shape of both matrices. 
 * @tparam C Column shape of both mtarices. 
 * @param x Left-hand side matrix. 
 * @param y Right-hand side matrix. 
 * @throw std::invalid_argument if sizes do not match. 
 */ 
template <typename LHS, typename RHS, int R, int C> 
inline void 
assign(Eigen::Matrix<LHS,R,C>& x, 
       const Eigen::Matrix<RHS,R,C>& y) { 
  if(x.rows() == 0 && x.cols() == 0) { 
    x = y; 
    return ; 
  } 
  stan::math::check_matching_dims("assign", 
                                            "x", x, 
                                            "y", y); 
  for (int i = 0; i < x.size(); ++i) 
    assign(x(i),y(i)); 
} 

[bgoodri]

Okay, it seems that

This instance will be called for arguments that are both Eigen::Matrix types and whose shapes match.

is a correct code comment for what the compiler interprets as a "match"; i.e. Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic> matches with another Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic> even if their rows and / or columns are different numbers. I think my hack of allowing resizing of 0x0 matrices also needs type promotion, but other than that I guess we only need the declaration change and I can PR a bunch more matrix decomposition stuff.

[bob-carpenter]

I'm still not convinced this is a good idea. I think it's
going to be confusing to users. I'm torn between wanting to
put in something useful versus holding off until we can
implement a proper solution. This is the same way I feel about
adding a single sparse matrix operation or very specific GLM
functions.

Is the plan to throw out size matches altogether in Stan's
assign and just allow Eigen-style assignment? Or to just allow
that if the lvalue is size 0? I'd be more comfortable with the
latter as it seems less error prone.

When you say "compiler", do you mean Stan's or C++'s? Eigen
doesn't do size matching. Stan's parser does not look at
size at all --- it just checks shape (e.g., matrix vs. vector,
1D array vs. scalar, etc.)

• Bob

On Mar 17, 2015, at 5:03 AM, bgoodri notifications@github.com wrote:

Okay, it seems that

This instance will be called for arguments that are both Eigen::Matrix types and whose shapes match.

is a correct code comment for what the compiler interprets as a "match"; i.e. Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic> matches with another Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic> even if their rows and / or columns are different numbers. I think my hack of allowing resizing of 0x0 matrices also needs type promotion, but other than that I guess we only need the declaration change and I can PR a bunch more matrix decomposition stuff.

—
Reply to this email directly or view it on GitHub.

[bgoodri]

I think having functions are able to return an array of matrices would be a good thing, and I am not seeing another way to accomplish that right now. Possibly we could make a push to allow simultaneous declaration and definitions, in which case things wouldn't really be "resized" but just given different assignments from the outset. But we would still need heterogenous sizes.

If we do resizing, it probably makes sense to only allow Eigen-style assignment if the lvalue has size 0 but it is I think more confusing to the user to have to write

matrix[0,0] Q_and_R[2];
Q_and_R <- QR(Sigma);

than to write

matrix Q_and_R[2];
Q_and_R <- QR(Sigma);

although

matrix Q_and_R[2] <- QR(Sigma);

would probably be better eventually. Although maybe when we go to C++, stanc could just use auto behind the scenes and then it could just be

Q_and_R <- QR(Sigma); // parses as auto Q_and_R = QR(Sigma);

The "compiler" thing that threw me was that the code comment was refering to Eigen's matching on Eigen::Dynamic whereas I was thinking in terms of algebra.

[bob-carpenter]

On Mar 17, 2015, at 10:28 AM, bgoodri notifications@github.com wrote:

I think having functions are able to return an array of matrices would be a good thing, and I am not seeing another way to accomplish that right now. Possibly we could make a push to allow simultaneous declaration and definitions,

That's on my near-term to-do list, but as you point out, it won't solve
the heterogeneous size problem.

in which case things wouldn't really be "resized" but just given different assignments from the outset. But we would still need heterogenous sizes.

Maybe we need an R list-like structure. The only
problem is how to declare such a thing and still allow
type inference.

If we do resizing, it probably makes sense to only allow Eigen-style assignment if the lvalue has size 0 but it is I think more confusing to the user to have to write

matrix[0,0] Q_and_R[2];
Q_R <- QR(Sigma);

than to write

matrix Q_and_R[2];
Q_and_R <- QR(Sigma);

Agreed. So that would involve a change in the language
parser and code generator. A simple way to deal with it
would be to take matrix as synonymous to "matrix[0,0]".

The user-defined function definitions also use bare
types, so there's already doc on what they look like.

We have to think in terms of generality here, so a user
could write

matrix a[3];

to get an array of matrix[0,0], each of which would then
be assignable to a different sized matrix.

We then have to think through what this does to things
like the num_elements and size functions.

although

matrix Q_and_R[2] <- QR(Sigma);

would probably be better eventually. Although maybe when we go to C++, stanc could just use auto behind the scenes and then it could just be

Q_and_R <- QR(Sigma); // parses as auto Q_and_R = QR(Sigma);

I don't like that idea because I want to preserve static typing,
as in C++, and not move to Python or R-style dynamic typing.

The "compiler" thing that threw me was that the code comment was refering to Eigen's matching on Eigen::Dynamic whereas I was thinking in terms of algebra.

Ah, right. You can use fixed sizes > 1, but they get implemented
as Dynamic at some small N, but tested statically for matching. I
decided not to do that in Stan because I thought it'd hugely complicate
argument passing to functions (built-ins, because I wasn't thinking
user-defined at the time).

• Bob

[bgoodri]

The closest thing to an R list is a std::tuple but then we would have to wait for C++11. I think if someone does

matrix Q_and_R[2];
print(size(Q_and_R)); // should still be 2
print(num_elements(Q_and_R)); // should still be 0
Q_and_R <- QR(Sigma);
print(num_elements(Q_and_R)); // should now be equal to rows(Sigma) * cols(Sigma)^3

[bob-carpenter]

Boost has already implemented tuples. They tend to
anticipate the new C++ work, since many of their package
authors are involved.

The problem I see is that if we have a list of size 100,
we'd need a size 100 typelist, which would just crush
a compiler.

Small tuples would be OK, though.

• Bob

On Mar 17, 2015, at 10:51 AM, bgoodri notifications@github.com wrote:

The closest thing to an R list is a std::tuple but then we would have to wait for C++11. I think if someone does

matrix Q_and_R[2];
print(size(Q_and_R)); // should still be 2
print(num_elements(Q_and_R)); // should still be 0
Q_and_R <- QR(Sigma);
print(num_elements(Q_and_R)); // should now be equal to rows(Sigma) * cols(Sigma)^3

—
Reply to this email directly or view it on GitHub.

[bgoodri]

I think tuples would be useful eventually (and hopefully the C++11 version is easier to compile), but they are more general than what is needed in this case where the heterogenous dynamic-sized matrices are "homogenous" from Eigen's perspective. If

matrix Q_and_R[2];

is equivalent to

matrix[0,0] Q_and_R[2];

it would work like user-defined functions and not disorient users much. I guess the danger is that someone acceidentally declares a 0x0 matrix and then assigns some other matrix to it like

matrix[K,K] foo; // where K = 0 due to a bug
foo <- bar; // where bar is not 0x0

But that might actually be the intended behavior and the user just meant for foo to me JxJ or something.

[bgoodri]

What is the right Stan syntax to do promotion inside here:

  if(x.rows() == 0 && x.cols() == 0) { 
    x = y;  // something like x = promote<scalar_type(x)>(y);
    return ; 
  } 

I need to implement this locally for a conference paper.

[bob-carpenter]

Is this in a user-defined function in Stan? It should
take care of whatever promotion you need automatically.

For instance, if you write this

real foo(real x, real y) {
return x + y;
}

and you pass in data (double) x and parameter (var) y,
the result will automatically be a var. And if x and y
are data and you assign to a local variable in the model
or transformed parameters block, you'll also get the
promotion for free.

• Bob

On Mar 19, 2015, at 6:13 AM, bgoodri notifications@github.com wrote:

What is the right Stan syntax to do promotion inside here:

if(x.rows() == 0 && x.cols() == 0) {
x = y; // something like x = promote<scalar_type(x)>(y);
return ;
}

I need to implement this locally for a conference paper.

—
Reply to this email directly or view it on GitHub.

[bgoodri]

No, I meant in assign.hpp to allow .stan programs with std::vectors of Eigen matrices with different sizes:

/** 
 * Copy the right-hand side's value to the left-hand side 
 * variable. 
 * 
 * The <code>assign()</code> function is overloaded.  This 
 * instance will be called for arguments that are both 
 * <code>Eigen::Matrix</code> types and whose shapes match.  The 
 * shapes are specified in the row and column template parameters. 
 * 
 * @tparam LHS Type of left-hand side matrix elements. 
 * @tparam RHS Type of right-hand side matrix elements. 
 * @tparam R Row shape of both matrices. 
 * @tparam C Column shape of both mtarices. 
 * @param x Left-hand side matrix. 
 * @param y Right-hand side matrix. 
 * @throw std::invalid_argument if sizes do not match. 
 */ 
template <typename LHS, typename RHS, int R, int C> 
inline void 
assign(Eigen::Matrix<LHS,R,C>& x, 
       const Eigen::Matrix<RHS,R,C>& y) { 
  if(x.rows() == 0 && x.cols() == 0) { 
    x = y; 
    return ; 
  } 
  stan::math::check_matching_dims("assign", 
                                            "x", x, 
                                            "y", y); 
  for (int i = 0; i < x.size(); ++i) 
    assign(x(i),y(i)); 
} 

[bob-carpenter]

Before going ahead with this, I'd like some feedback from
everyone else that they think it's a good idea to allow:

matrix[2,3] y;
matrix[0,0] x;
x <- y;

or

matrix[2,3] y;
matrix x;
x <- y;

I'm less happy with the former, but if we allow the latter,
we're pretty much going to have to allow the former.

• Bob

On Mar 19, 2015, at 10:35 AM, bgoodri notifications@github.com wrote:

No, I meant in assign.hpp to allow .stan programs with std::vectors of Eigen matrices with different sizes:

/**

• Copy the right-hand side's value to the left-hand side
• variable.
• The assign() function is overloaded. This
• instance will be called for arguments that are both
• Eigen::Matrix types and whose shapes match. The
• shapes are specified in the row and column template parameters.
• @tparam LHS Type of left-hand side matrix elements.
• @tparam RHS Type of right-hand side matrix elements.
• @tparam R Row shape of both matrices.
• @tparam C Column shape of both mtarices.
• @param x Left-hand side matrix.
• @param y Right-hand side matrix.
• @throw std::invalid_argument if sizes do not match.
  */
  template <typename LHS, typename RHS, int R, int C>
  inline void
  assign(Eigen::Matrix<LHS,R,C>& x,
  const Eigen::Matrix<RHS,R,C>& y) {
  if(x.rows() == 0 && x.cols() == 0) {
  x = y;
  return ;
  }
  stan::math::check_matching_dims("assign",
  "x", x,
  "y", y);
  for (int i = 0; i < x.size(); ++i)
  assign(x(i),y(i));
  }

—
Reply to this email directly or view it on GitHub.

[bob-carpenter]

Do you need to do anything other than
change check_matching_dims() to allow size zero
matrices or vectors on the LHS to be reassigned?

This is all going to be a bit confusing because we'll
allow arrays of vectors/matrices of different sizes,
but not arrays of arrays of different sizes.

• Bob

On Mar 19, 2015, at 10:35 AM, bgoodri notifications@github.com wrote:

No, I meant in assign.hpp to allow .stan programs with std::vectors of Eigen matrices with different sizes:

/**

• Copy the right-hand side's value to the left-hand side
• variable.
• The assign() function is overloaded. This
• instance will be called for arguments that are both
• Eigen::Matrix types and whose shapes match. The
• shapes are specified in the row and column template parameters.
• @tparam LHS Type of left-hand side matrix elements.
• @tparam RHS Type of right-hand side matrix elements.
• @tparam R Row shape of both matrices.
• @tparam C Column shape of both mtarices.
• @param x Left-hand side matrix.
• @param y Right-hand side matrix.
• @throw std::invalid_argument if sizes do not match.
  */
  template <typename LHS, typename RHS, int R, int C>
  inline void
  assign(Eigen::Matrix<LHS,R,C>& x,
  const Eigen::Matrix<RHS,R,C>& y) {
  if(x.rows() == 0 && x.cols() == 0) {
  x = y;
  return ;
  }
  stan::math::check_matching_dims("assign",
  "x", x,
  "y", y);
  for (int i = 0; i < x.size(); ++i)
  assign(x(i),y(i));
  }

—
Reply to this email directly or view it on GitHub.

[bgoodri]

Maybe. But check_matching_dims is called in several places besides
assign.hpp:

goodrich@T540p:/opt/stan$ git grep -F -n "check_matching_dims" src/stan |
cat
src/stan/math/fwd/mat/fun/columns_dot_product.hpp:7:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/fwd/mat/fun/columns_dot_product.hpp:20:
stan::math::check_matching_dims("columns_dot_product",
src/stan/math/fwd/mat/fun/columns_dot_product.hpp:37:
stan::math::check_matching_dims("columns_dot_product",
src/stan/math/fwd/mat/fun/columns_dot_product.hpp:54:
stan::math::check_matching_dims("columns_dot_product",
src/stan/math/fwd/mat/fun/rows_dot_product.hpp:7:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/fwd/mat/fun/rows_dot_product.hpp:21:
stan::math::check_matching_dims("rows_dot_product",
src/stan/math/fwd/mat/fun/rows_dot_product.hpp:38:
stan::math::check_matching_dims("rows_dot_product",
src/stan/math/fwd/mat/fun/rows_dot_product.hpp:55:
stan::math::check_matching_dims("rows_dot_product",
src/stan/math/fwd/mat/fun/to_fvar.hpp:6:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/fwd/mat/fun/to_fvar.hpp:51:
stan::math::check_matching_dims("to_fvar",
src/stan/math/prim/mat.hpp:18:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/prim/mat/err/check_matching_dims.hpp:37: inline bool
check_matching_dims(const char* function,
src/stan/math/prim/mat/fun/add.hpp:6:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/prim/mat/fun/add.hpp:29:
stan::math::check_matching_dims("add",
src/stan/math/prim/mat/fun/assign.hpp:11:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/prim/mat/fun/assign.hpp:98:
stan::math::check_matching_dims("assign",
src/stan/math/prim/mat/fun/elt_divide.hpp:6:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/prim/mat/fun/elt_divide.hpp:26:
stan::math::check_matching_dims("elt_divide",
src/stan/math/prim/mat/fun/elt_multiply.hpp:6:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/prim/mat/fun/elt_multiply.hpp:27:
stan::math::check_matching_dims("elt_multiply",
src/stan/math/prim/mat/fun/subtract.hpp:6:#include
<stan/math/prim/mat/err/check_matching_dims.hpp>
src/stan/math/prim/mat/fun/subtract.hpp:29:
stan::math::check_matching_dims("subtract",

and I didn't want to create problems if 0x0 matrices are passed to any of
those.

On Wed, Mar 18, 2015 at 8:00 PM, Bob Carpenter notifications@github.com
wrote:

Do you need to do anything other than
change check_matching_dims() to allow size zero
matrices or vectors on the LHS to be reassigned?

This is all going to be a bit confusing because we'll
allow arrays of vectors/matrices of different sizes,
but not arrays of arrays of different sizes.

• Bob

On Mar 19, 2015, at 10:35 AM, bgoodri notifications@github.com wrote:

No, I meant in assign.hpp to allow .stan programs with std::vectors of
Eigen matrices with different sizes:

/**

• Copy the right-hand side's value to the left-hand side
• variable.
  *
• The assign() function is overloaded. This
• instance will be called for arguments that are both
• Eigen::Matrix types and whose shapes match. The
• shapes are specified in the row and column template parameters.
  *
• @tparam LHS Type of left-hand side matrix elements.
• @tparam RHS Type of right-hand side matrix elements.
• @tparam R Row shape of both matrices.
• @tparam C Column shape of both mtarices.
• @param x Left-hand side matrix.
• @param y Right-hand side matrix.
• @throw std::invalid_argument if sizes do not match.
  */
  template <typename LHS, typename RHS, int R, int C>
  inline void
  assign(Eigen::Matrix<LHS,R,C>& x,
  const Eigen::Matrix<RHS,R,C>& y) {
  if(x.rows() == 0 && x.cols() == 0) {
  x = y;
  return ;
  }
  stan::math::check_matching_dims("assign",
  "x", x,
  "y", y);
  for (int i = 0; i < x.size(); ++i)
  assign(x(i),y(i));
  }

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bgoodri]

I agree that it seems hard to disallow the former.

But a matrix[0,0] is currently not useful for anything substantive; it just
means that you don't have to change your code when you leave a variable out
of a model. So, making matrix[0,0] tantamount to "unknown rows and columns"
is a small step forward in the usefulness department, it just looks worse
than an undimensioned matrix declaration so I doubt it would often get
written intentionally. It may rarely get written unintentionally but that
is a self-inflicted bug by the user that would usually get exposed later in
the Stan code.

On Wed, Mar 18, 2015 at 7:56 PM, Bob Carpenter notifications@github.com
wrote:

Before going ahead with this, I'd like some feedback from
everyone else that they think it's a good idea to allow:

matrix[2,3] y;
matrix[0,0] x;
x <- y;

or

matrix[2,3] y;
matrix x;
x <- y;

I'm less happy with the former, but if we allow the latter,
we're pretty much going to have to allow the former.

• Bob

On Mar 19, 2015, at 10:35 AM, bgoodri notifications@github.com wrote:

No, I meant in assign.hpp to allow .stan programs with std::vectors of
Eigen matrices with different sizes:

/**

• Copy the right-hand side's value to the left-hand side
• variable.
  *
• The assign() function is overloaded. This
• instance will be called for arguments that are both
• Eigen::Matrix types and whose shapes match. The
• shapes are specified in the row and column template parameters.
  *
• @tparam LHS Type of left-hand side matrix elements.
• @tparam RHS Type of right-hand side matrix elements.
• @tparam R Row shape of both matrices.
• @tparam C Column shape of both mtarices.
• @param x Left-hand side matrix.
• @param y Right-hand side matrix.
• @throw std::invalid_argument if sizes do not match.
  */
  template <typename LHS, typename RHS, int R, int C>
  inline void
  assign(Eigen::Matrix<LHS,R,C>& x,
  const Eigen::Matrix<RHS,R,C>& y) {
  if(x.rows() == 0 && x.cols() == 0) {
  x = y;
  return ;
  }
  stan::math::check_matching_dims("assign",
  "x", x,
  "y", y);
  for (int i = 0; i < x.size(); ++i)
  assign(x(i),y(i));
  }

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[syclik]

Sorry for taking so long to catch up.

First off, I disagree with allowing matrix[0,0] stand for a resizable matrix. Hopefully Stan is now compatible with 0 x 0 matrices, 0-length vectors, and 0-length row vectors. Early on, I found myself wanting to draw from priors and would intentionally write models with size 0 to allow for this. Perhaps the behavior to use this functionality within the Stan language wouldn't be affected, but I think this might allow users (me) to write models that behaved differently than intended.

Regarding matrix without sizes representing a resizable matrix, it does make some sense, but I'm afraid it'll make the user experience a little worse. The problem is when people forget to size matrices, we'll lose the power to write better error messages due to the strict typing at that point. Yes, I know, we don't proactively stop now, but the code is almost ready to do so.

Rather than just saying no, my suggestion would be to create a new type, something called dynamic_matrix. (I was inspired by Eigen's naming... and this is just a suggestion, we can call it whatever makes sense.) I'm really trying not to get ahead of myself and thinking about implementation because the C++ types would be the same (Eigen::Matrix), but we would want to do different checking. Perhaps we would need to have another set of assign calls or something like that? Not sure.

In short: I don't like matrix[0,0] or just matrix, but would prefer dynamic_matrix.

[bgoodri]

This wouldn't limit the ability to use a size zero matrix as a placeholder
in the model. If you don't resize it, then it stays size zero.

I'm not sure what the code is almost ready to do with regards to error
messages. I think what would happen is whatever would happen if someone
declared a size zero matrix and then did something buggy:

matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
X <- rep_matrix(0, 2, 2); // resized to 2x2
print(X[1,2]); // prints 0 correctly

I don't think adding the Stan keyword dynamic would hurt anything but how
it help if the construct gets parsed to an
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> in C++? Are you thinking of
subclassing from Eigen::Matrix?

On Fri, Mar 20, 2015 at 10:46 AM, Daniel Lee notifications@github.com
wrote:

Sorry for taking so long to catch up.

First off, I disagree with allowing matrix[0,0] stand for a resizable
matrix. Hopefully Stan is now compatible with 0 x 0 matrices, 0-length
vectors, and 0-length row vectors. Early on, I found myself wanting to draw
from priors and would intentionally write models with size 0 to allow for
this. Perhaps the behavior to use this functionality within the Stan
language wouldn't be affected, but I think this might allow users (me) to
write models that behaved differently than intended.

Regarding matrix without sizes representing a resizable matrix, it does
make some sense, but I'm afraid it'll make the user experience a little
worse. The problem is when people forget to size matrices, we'll lose the
power to write better error messages due to the strict typing at that
point. Yes, I know, we don't proactively stop now, but the code is almost
ready to do so.

Rather than just saying no, my suggestion would be to create a new type,
something called dynamic_matrix. (I was inspired by Eigen's naming... and
this is just a suggestion, we can call it whatever makes sense.) I'm really
trying not to get ahead of myself and thinking about implementation because
the C++ types would be the same (Eigen::Matrix), but we would want to do
different checking. Perhaps we would need to have another set of assign
calls or something like that? Not sure.

In short: I don't like matrix[0,0] or just matrix, but would prefer
dynamic_matrix.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[betanalpha]

I agree that if a type is meant to be resizable
we should make that explicit. I think most
people assume that their matrix sizes are
constant, and that’s a good thing.

On Mar 20, 2015, at 3:25 PM, bgoodri notifications@github.com wrote:

This wouldn't limit the ability to use a size zero matrix as a placeholder
in the model. If you don't resize it, then it stays size zero.

I'm not sure what the code is almost ready to do with regards to error
messages. I think what would happen is whatever would happen if someone
declared a size zero matrix and then did something buggy:

matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
X <- rep_matrix(0, 2, 2); // resized to 2x2
print(X[1,2]); // prints 0 correctly

I don't think adding the Stan keyword dynamic would hurt anything but how
it help if the construct gets parsed to an
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> in C++? Are you thinking of
subclassing from Eigen::Matrix?

On Fri, Mar 20, 2015 at 10:46 AM, Daniel Lee notifications@github.com
wrote:

Sorry for taking so long to catch up.

First off, I disagree with allowing matrix[0,0] stand for a resizable
matrix. Hopefully Stan is now compatible with 0 x 0 matrices, 0-length
vectors, and 0-length row vectors. Early on, I found myself wanting to draw
from priors and would intentionally write models with size 0 to allow for
this. Perhaps the behavior to use this functionality within the Stan
language wouldn't be affected, but I think this might allow users (me) to
write models that behaved differently than intended.

Regarding matrix without sizes representing a resizable matrix, it does
make some sense, but I'm afraid it'll make the user experience a little
worse. The problem is when people forget to size matrices, we'll lose the
power to write better error messages due to the strict typing at that
point. Yes, I know, we don't proactively stop now, but the code is almost
ready to do so.

Rather than just saying no, my suggestion would be to create a new type,
something called dynamic_matrix. (I was inspired by Eigen's naming... and
this is just a suggestion, we can call it whatever makes sense.) I'm really
trying not to get ahead of myself and thinking about implementation because
the C++ types would be the same (Eigen::Matrix), but we would want to do
different checking. Perhaps we would need to have another set of assign
calls or something like that? Not sure.

In short: I don't like matrix[0,0] or just matrix, but would prefer
dynamic_matrix.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

[bgoodri]

Do you mean, just for users to remind themselves that it is resizeable?

resizeable matrix X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

That's fine by me, but how would it be different on the C++ as

matrix[0,0] X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

On Fri, Mar 20, 2015 at 11:49 AM, Michael Betancourt <
notifications@github.com> wrote:

I agree that if a type is meant to be resizable
we should make that explicit. I think most
people assume that their matrix sizes are
constant, and that’s a good thing.

On Mar 20, 2015, at 3:25 PM, bgoodri notifications@github.com wrote:

This wouldn't limit the ability to use a size zero matrix as a
placeholder
in the model. If you don't resize it, then it stays size zero.

I'm not sure what the code is almost ready to do with regards to error
messages. I think what would happen is whatever would happen if someone
declared a size zero matrix and then did something buggy:

matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
X <- rep_matrix(0, 2, 2); // resized to 2x2
print(X[1,2]); // prints 0 correctly

I don't think adding the Stan keyword dynamic would hurt anything but how
it help if the construct gets parsed to an
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> in C++? Are you thinking
of
subclassing from Eigen::Matrix?

On Fri, Mar 20, 2015 at 10:46 AM, Daniel Lee notifications@github.com
wrote:

Sorry for taking so long to catch up.

First off, I disagree with allowing matrix[0,0] stand for a resizable
matrix. Hopefully Stan is now compatible with 0 x 0 matrices, 0-length
vectors, and 0-length row vectors. Early on, I found myself wanting to
draw
from priors and would intentionally write models with size 0 to allow
for
this. Perhaps the behavior to use this functionality within the Stan
language wouldn't be affected, but I think this might allow users (me)
to
write models that behaved differently than intended.

Regarding matrix without sizes representing a resizable matrix, it does
make some sense, but I'm afraid it'll make the user experience a little
worse. The problem is when people forget to size matrices, we'll lose
the
power to write better error messages due to the strict typing at that
point. Yes, I know, we don't proactively stop now, but the code is
almost
ready to do so.

Rather than just saying no, my suggestion would be to create a new
type,
something called dynamic_matrix. (I was inspired by Eigen's naming...
and
this is just a suggestion, we can call it whatever makes sense.) I'm
really
trying not to get ahead of myself and thinking about implementation
because
the C++ types would be the same (Eigen::Matrix), but we would want to
do
different checking. Perhaps we would need to have another set of assign
calls or something like that? Not sure.

In short: I don't like matrix[0,0] or just matrix, but would prefer
dynamic_matrix.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[betanalpha]

Yeah, I’m just thinking about the users — the
underlying code can be exactly the same.

On Mar 20, 2015, at 3:58 PM, bgoodri notifications@github.com wrote:

Do you mean, just for users to remind themselves that it is resizeable?

resizeable matrix X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

That's fine by me, but how would it be different on the C++ as

matrix[0,0] X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

On Fri, Mar 20, 2015 at 11:49 AM, Michael Betancourt <
notifications@github.com> wrote:

I agree that if a type is meant to be resizable
we should make that explicit. I think most
people assume that their matrix sizes are
constant, and that’s a good thing.

On Mar 20, 2015, at 3:25 PM, bgoodri notifications@github.com wrote:

This wouldn't limit the ability to use a size zero matrix as a
placeholder
in the model. If you don't resize it, then it stays size zero.

I'm not sure what the code is almost ready to do with regards to error
messages. I think what would happen is whatever would happen if someone
declared a size zero matrix and then did something buggy:

matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
X <- rep_matrix(0, 2, 2); // resized to 2x2
print(X[1,2]); // prints 0 correctly

I don't think adding the Stan keyword dynamic would hurt anything but how
it help if the construct gets parsed to an
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> in C++? Are you thinking
of
subclassing from Eigen::Matrix?

On Fri, Mar 20, 2015 at 10:46 AM, Daniel Lee notifications@github.com
wrote:

Sorry for taking so long to catch up.

First off, I disagree with allowing matrix[0,0] stand for a resizable
matrix. Hopefully Stan is now compatible with 0 x 0 matrices, 0-length
vectors, and 0-length row vectors. Early on, I found myself wanting to
draw
from priors and would intentionally write models with size 0 to allow
for
this. Perhaps the behavior to use this functionality within the Stan
language wouldn't be affected, but I think this might allow users (me)
to
write models that behaved differently than intended.

Regarding matrix without sizes representing a resizable matrix, it does
make some sense, but I'm afraid it'll make the user experience a little
worse. The problem is when people forget to size matrices, we'll lose
the
power to write better error messages due to the strict typing at that
point. Yes, I know, we don't proactively stop now, but the code is
almost
ready to do so.

Rather than just saying no, my suggestion would be to create a new
type,
something called dynamic_matrix. (I was inspired by Eigen's naming...
and
this is just a suggestion, we can call it whatever makes sense.) I'm
really
trying not to get ahead of myself and thinking about implementation
because
the C++ types would be the same (Eigen::Matrix), but we would want to
do
different checking. Perhaps we would need to have another set of assign
calls or something like that? Not sure.

In short: I don't like matrix[0,0] or just matrix, but would prefer
dynamic_matrix.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

[bgoodri]

If Bob wants to add a dynamic or resizeable keyword, go ahead. One concern
is that if someone declares a resizeable matrix, they might think that they
can resize it multiple times, when really it is more of a "deferred size"
thing: it is temporarily 0x0 and gets changed to a fixed size that is not
0x0.

On Fri, Mar 20, 2015 at 12:00 PM, Michael Betancourt <
notifications@github.com> wrote:

Yeah, I’m just thinking about the users — the
underlying code can be exactly the same.

On Mar 20, 2015, at 3:58 PM, bgoodri notifications@github.com wrote:

Do you mean, just for users to remind themselves that it is resizeable?

resizeable matrix X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

That's fine by me, but how would it be different on the C++ as

matrix[0,0] X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

On Fri, Mar 20, 2015 at 11:49 AM, Michael Betancourt <
notifications@github.com> wrote:

I agree that if a type is meant to be resizable
we should make that explicit. I think most
people assume that their matrix sizes are
constant, and that’s a good thing.

On Mar 20, 2015, at 3:25 PM, bgoodri notifications@github.com wrote:

This wouldn't limit the ability to use a size zero matrix as a
placeholder
in the model. If you don't resize it, then it stays size zero.

I'm not sure what the code is almost ready to do with regards to
error
messages. I think what would happen is whatever would happen if
someone
declared a size zero matrix and then did something buggy:

matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
X <- rep_matrix(0, 2, 2); // resized to 2x2
print(X[1,2]); // prints 0 correctly

I don't think adding the Stan keyword dynamic would hurt anything
but how
it help if the construct gets parsed to an
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> in C++? Are you
thinking
of
subclassing from Eigen::Matrix?

On Fri, Mar 20, 2015 at 10:46 AM, Daniel Lee <
notifications@github.com>
wrote:

Sorry for taking so long to catch up.

First off, I disagree with allowing matrix[0,0] stand for a
resizable
matrix. Hopefully Stan is now compatible with 0 x 0 matrices,
0-length
vectors, and 0-length row vectors. Early on, I found myself
wanting to
draw
from priors and would intentionally write models with size 0 to
allow
for
this. Perhaps the behavior to use this functionality within the
Stan
language wouldn't be affected, but I think this might allow users
(me)
to
write models that behaved differently than intended.

Regarding matrix without sizes representing a resizable matrix, it
does
make some sense, but I'm afraid it'll make the user experience a
little
worse. The problem is when people forget to size matrices, we'll
lose
the
power to write better error messages due to the strict typing at
that
point. Yes, I know, we don't proactively stop now, but the code is
almost
ready to do so.

Rather than just saying no, my suggestion would be to create a new
type,
something called dynamic_matrix. (I was inspired by Eigen's
naming...
and
this is just a suggestion, we can call it whatever makes sense.)
I'm
really
trying not to get ahead of myself and thinking about implementation
because
the C++ types would be the same (Eigen::Matrix), but we would want
to
do
different checking. Perhaps we would need to have another set of
assign
calls or something like that? Not sure.

In short: I don't like matrix[0,0] or just matrix, but would prefer
dynamic_matrix.

—
Reply to this email directly or view it on GitHub
<
https://github.com/stan-dev/stan/issues/1369#issuecomment-84034932>.

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bob-carpenter]

Adding a no-op keyword is easy, but I don't think that's
a good idea.

Right now, we allow size-zero vectors and matrices, they print
just fine, and can be used as data to provide a no-op in
appropriate situations (e.g., no data, or no predictions).

So Ben's right that if we hack in a solution by allowing
matrix[0,0](and presumably vector[0] and row_vector[0] for
consistency) to be assigned to anything, then

• they would not be reassignable after the first non-zero assign, and
• a matrix declared as matrix[0,0] would still be reassigable

Anything else is going to take more work, like adding
a truly resizable matrix data type. I'm not even sure how
to do that with the code generation we have other than creating
a whole new type. And even then, I'm not sure.

I was thinking we'd have ragged arrays with declared array sizes.
Something like

int Ms[K];
int Ns[K];
matrix[Ms,Ns] a[K];

where a[k] is of size Ms[k] x Ns[k]. That's consistent with how
we do things now --- all sizes fixed at declaration time other than
function arguments and returns, which are fixed on call and
return rather than when declared.

• Bob

On Mar 21, 2015, at 3:06 AM, bgoodri notifications@github.com wrote:

If Bob wants to add a dynamic or resizeable keyword, go ahead. One concern
is that if someone declares a resizeable matrix, they might think that they
can resize it multiple times, when really it is more of a "deferred size"
thing: it is temporarily 0x0 and gets changed to a fixed size that is not
0x0.

On Fri, Mar 20, 2015 at 12:00 PM, Michael Betancourt <
notifications@github.com> wrote:

Yeah, I’m just thinking about the users — the
underlying code can be exactly the same.

On Mar 20, 2015, at 3:58 PM, bgoodri notifications@github.com wrote:

Do you mean, just for users to remind themselves that it is resizeable?

resizeable matrix X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

That's fine by me, but how would it be different on the C++ as

matrix[0,0] X;
X <- rep_matrix(0, 2, 2);
print(X[1,2]);

On Fri, Mar 20, 2015 at 11:49 AM, Michael Betancourt <
notifications@github.com> wrote:

I agree that if a type is meant to be resizable
we should make that explicit. I think most
people assume that their matrix sizes are
constant, and that’s a good thing.

On Mar 20, 2015, at 3:25 PM, bgoodri notifications@github.com wrote:

This wouldn't limit the ability to use a size zero matrix as a
placeholder
in the model. If you don't resize it, then it stays size zero.

I'm not sure what the code is almost ready to do with regards to
error
messages. I think what would happen is whatever would happen if
someone
declared a size zero matrix and then did something buggy:

matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
print(X[1,2]); // error

matrix X; // equivalent to matrix[0,0] X;
X <- rep_matrix(0, 2, 2); // resized to 2x2
print(X[1,2]); // prints 0 correctly

I don't think adding the Stan keyword dynamic would hurt anything
but how
it help if the construct gets parsed to an
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> in C++? Are you
thinking
of
subclassing from Eigen::Matrix?

On Fri, Mar 20, 2015 at 10:46 AM, Daniel Lee <
notifications@github.com>
wrote:

Sorry for taking so long to catch up.

First off, I disagree with allowing matrix[0,0] stand for a
resizable
matrix. Hopefully Stan is now compatible with 0 x 0 matrices,
0-length
vectors, and 0-length row vectors. Early on, I found myself
wanting to
draw
from priors and would intentionally write models with size 0 to
allow
for
this. Perhaps the behavior to use this functionality within the
Stan
language wouldn't be affected, but I think this might allow users
(me)
to
write models that behaved differently than intended.

Regarding matrix without sizes representing a resizable matrix, it
does
make some sense, but I'm afraid it'll make the user experience a
little
worse. The problem is when people forget to size matrices, we'll
lose
the
power to write better error messages due to the strict typing at
that
point. Yes, I know, we don't proactively stop now, but the code is
almost
ready to do so.

Rather than just saying no, my suggestion would be to create a new
type,
something called dynamic_matrix. (I was inspired by Eigen's
naming...
and
this is just a suggestion, we can call it whatever makes sense.)
I'm
really
trying not to get ahead of myself and thinking about implementation
because
the C++ types would be the same (Eigen::Matrix), but we would want
to
do
different checking. Perhaps we would need to have another set of
assign
calls or something like that? Not sure.

In short: I don't like matrix[0,0] or just matrix, but would prefer
dynamic_matrix.

—
Reply to this email directly or view it on GitHub
<
https://github.com/stan-dev/stan/issues/1369#issuecomment-84034932>.

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

[syclik]

I think Bob's suggestion is what makes the most sense. I'm guessing will take a long time to figure out, so, let's talk about what we can do in the short term.

@bgoodri, I still think allowing matrix[0,0] to mean something other than a zero-row, zero-column matrix is a bit dangerous and could lead to model building errors that aren't clear. I also think having a naked matrix call is also just asking for trouble, but it's better in my book.

---

I know, I didn't come up with a decent counter-example, but perhaps something like:

data {
   int N;
   int M;
   matrix[N,N] X1;
   matrix[M,M] X2;
}
transformed data {
   matrix[M,M] L_X1;
   matrix[M,M] L_X2;
   L_X1 <- cholesky_decompose(X1);
   L_X2 <- cholesky_decompose(X2);
}
...

That'll work if M is 0. Seems odd and error-prone to me.

[bgoodri]

If we implement true ragged arrays, that would be better because it would
work in data and parameters as well. But we've been talking about doing
that since 2011. And only Bob can do something like that with the language.
Whereas it takes 4 lines of already-written code to assign to a 0x0 matrix
to overcome a restriction we imposed on the language for a conference paper.

On Wed, Mar 25, 2015 at 3:13 PM, Daniel Lee notifications@github.com
wrote:

I think Bob's suggestion is what makes the most sense. I'm guessing will
take a long time to figure out, so, let's talk about what we can do in the
short term.

@bgoodri https://github.com/bgoodri, I still think allowing matrix[0,0]
to mean something other than a zero-row, zero-column matrix is a bit
dangerous and could lead to model building errors that aren't clear. I also
think having a naked matrix call is also just asking for trouble, but

it's better in my book.

I know, I didn't come up with a decent counter-example, but perhaps
something like:

data {
int N;
int M;
matrix[N,N] X1;
matrix[M,M] X2;
}
transformed data {
matrix[M,M] L_X1;
matrix[M,M] L_X2;
L_X1 <- cholesky_decompose(X1);
L_X2 <- cholesky_decompose(X2);
}
...

That'll work if M is 0. Seems odd and error-prone to me.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bob-carpenter]

This project has to be a compromise. I don't like this
proposal because it messes up the clean conception of sizing
that we have everywhere. But I agree that we've been wanting
to do this forever. If I could ever stop treading water, I
could get lots of these features in. Maybe I'll try not responding
to the mailing list again and see if anyone picks up the slack.
Daniel's been volunteering.

Don't forget doc and testing. It's may be 4 lines of
code to change the assign function, but then there's testing
and mods to the doc everywhere we talk about assignment.

By the way, is this only for matrices, or also for vectors
and row vectors?

I think it'll actually require a fair amount of doc and testing.
Can we do this reassignment in the transformed data block or just
with local variables? If we can do it in the former, then
what happens to the validation tests for sizes that run at the
end of the blocks --- those will also have to change to match with
tests to go along? I don't think we can do it in transformed parameters,
because that could change the number of parameters on a per-iteration
basis. So then all the discussion of matrix assignment in the doc
needs to be qualified.

Maybe we can come up with an intermediate solution, such as having
a matrix declared without any size restrictions, as in

matrix A;

And then only use this as a local variable. Or maybe in transformed
data.

Then when using that on the LHS of an array, we'd not check sizes.
That'd require a lot more code and it'd change explicitly the
sort of implicit change that Ben's arguing for. That is, we'd now
have a contrast between sized and unsized types, allowing assignment
and presumably reassignment to unsized types. That's even more code
and doc, but I think it's a much more coherent solution than modifying
the behavior of size-0 matrices.

• Bob

On Mar 26, 2015, at 8:59 AM, bgoodri notifications@github.com wrote:

If we implement true ragged arrays, that would be better because it would
work in data and parameters as well. But we've been talking about doing
that since 2011. And only Bob can do something like that with the language.
Whereas it takes 4 lines of already-written code to assign to a 0x0 matrix
to overcome a restriction we imposed on the language for a conference paper.

On Wed, Mar 25, 2015 at 3:13 PM, Daniel Lee notifications@github.com
wrote:

I think Bob's suggestion is what makes the most sense. I'm guessing will
take a long time to figure out, so, let's talk about what we can do in the
short term.

@bgoodri https://github.com/bgoodri, I still think allowing matrix[0,0]
to mean something other than a zero-row, zero-column matrix is a bit
dangerous and could lead to model building errors that aren't clear. I also
think having a naked matrix call is also just asking for trouble, but

it's better in my book.

I know, I didn't come up with a decent counter-example, but perhaps
something like:

data {
int N;
int M;
matrix[N,N] X1;
matrix[M,M] X2;
}
transformed data {
matrix[M,M] L_X1;
matrix[M,M] L_X2;
L_X1 <- cholesky_decompose(X1);
L_X2 <- cholesky_decompose(X2);
}
...

That'll work if M is 0. Seems odd and error-prone to me.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

[bgoodri]

If we are going to support genuine ragged arrays like

matrix[rows,cols] X[elements];

where rows and cols are vectors, then I think we should allow that anywhere.

If it is just the original hack where a 0x0 matrix or 0x1 vector can get
reassigned to any genuine size, then it definitely can't be data or
parameters and if it can only be model or generated quantities, then that
would still be useful.

More generally, I am not a fan of limiting language features that lots of
people would find useful in order to catch bugs in user code. If we can
catch bugs like that essentially costlessly, then okay. But, such bugs will
probably manifest themselves elsewhere and if not, too bad for that
particular user, as opposed to bad for all users. So, in this case, are you
saying that

transformed parameters {
cov_matrix[K,K] X;
X <- crossprod(Y);
}

currently checks at the end of the block that X is, in fact, KxK. If so, I
agree it would be hard to check if X were declared 0x0 and resized, but I
would favor dropping that check. But we could still check that it is
positive definite without regard to the size, right?

On Wed, Mar 25, 2015 at 7:05 PM, Bob Carpenter notifications@github.com
wrote:

This project has to be a compromise. I don't like this
proposal because it messes up the clean conception of sizing
that we have everywhere. But I agree that we've been wanting
to do this forever. If I could ever stop treading water, I
could get lots of these features in. Maybe I'll try not responding
to the mailing list again and see if anyone picks up the slack.
Daniel's been volunteering.

Don't forget doc and testing. It's may be 4 lines of
code to change the assign function, but then there's testing
and mods to the doc everywhere we talk about assignment.

By the way, is this only for matrices, or also for vectors
and row vectors?

I think it'll actually require a fair amount of doc and testing.
Can we do this reassignment in the transformed data block or just
with local variables? If we can do it in the former, then
what happens to the validation tests for sizes that run at the
end of the blocks --- those will also have to change to match with
tests to go along? I don't think we can do it in transformed parameters,
because that could change the number of parameters on a per-iteration
basis. So then all the discussion of matrix assignment in the doc
needs to be qualified.

Maybe we can come up with an intermediate solution, such as having
a matrix declared without any size restrictions, as in

matrix A;

And then only use this as a local variable. Or maybe in transformed
data.

Then when using that on the LHS of an array, we'd not check sizes.
That'd require a lot more code and it'd change explicitly the
sort of implicit change that Ben's arguing for. That is, we'd now
have a contrast between sized and unsized types, allowing assignment
and presumably reassignment to unsized types. That's even more code
and doc, but I think it's a much more coherent solution than modifying
the behavior of size-0 matrices.

• Bob

On Mar 26, 2015, at 8:59 AM, bgoodri notifications@github.com wrote:

If we implement true ragged arrays, that would be better because it would
work in data and parameters as well. But we've been talking about doing
that since 2011. And only Bob can do something like that with the
language.
Whereas it takes 4 lines of already-written code to assign to a 0x0
matrix
to overcome a restriction we imposed on the language for a conference
paper.

On Wed, Mar 25, 2015 at 3:13 PM, Daniel Lee notifications@github.com
wrote:

I think Bob's suggestion is what makes the most sense. I'm guessing
will
take a long time to figure out, so, let's talk about what we can do in
the
short term.

@bgoodri https://github.com/bgoodri, I still think allowing
matrix[0,0]
to mean something other than a zero-row, zero-column matrix is a bit
dangerous and could lead to model building errors that aren't clear. I
also
think having a naked matrix call is also just asking for trouble, but

it's better in my book.

I know, I didn't come up with a decent counter-example, but perhaps
something like:

data {
int N;
int M;
matrix[N,N] X1;
matrix[M,M] X2;
}
transformed data {
matrix[M,M] L_X1;
matrix[M,M] L_X2;
L_X1 <- cholesky_decompose(X1);
L_X2 <- cholesky_decompose(X2);
}
...

That'll work if M is 0. Seems odd and error-prone to me.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bob-carpenter]

I agree about allowing declared ragged arrays anywhere. And I
fully admit I don't have this on my to-do list for any time soon.

I strongly disagree about turning off size checks at the ends of
blocks. And in general about favoring functionality over error
checking.

Maybe we can discuss at the next Stan meeting?

On Mar 26, 2015, at 10:27 AM, bgoodri notifications@github.com wrote:

If we are going to support genuine ragged arrays like

matrix[rows,cols] X[elements];

where rows and cols are vectors, then I think we should allow that anywhere.

If it is just the original hack where a 0x0 matrix or 0x1 vector can get
reassigned to any genuine size, then it definitely can't be data or
parameters and if it can only be model or generated quantities, then that
would still be useful.

More generally, I am not a fan of limiting language features that lots of
people would find useful in order to catch bugs in user code. If we can
catch bugs like that essentially costlessly, then okay. But, such bugs will
probably manifest themselves elsewhere and if not, too bad for that
particular user, as opposed to bad for all users. So, in this case, are you
saying that

transformed parameters {
cov_matrix[K,K] X;
X <- crossprod(Y);
}

currently checks at the end of the block that X is, in fact, KxK. If so, I
agree it would be hard to check if X were declared 0x0 and resized, but I
would favor dropping that check. But we could still check that it is
positive definite without regard to the size, right?

On Wed, Mar 25, 2015 at 7:05 PM, Bob Carpenter notifications@github.com
wrote:

This project has to be a compromise. I don't like this
proposal because it messes up the clean conception of sizing
that we have everywhere. But I agree that we've been wanting
to do this forever. If I could ever stop treading water, I
could get lots of these features in. Maybe I'll try not responding
to the mailing list again and see if anyone picks up the slack.
Daniel's been volunteering.

Don't forget doc and testing. It's may be 4 lines of
code to change the assign function, but then there's testing
and mods to the doc everywhere we talk about assignment.

By the way, is this only for matrices, or also for vectors
and row vectors?

I think it'll actually require a fair amount of doc and testing.
Can we do this reassignment in the transformed data block or just
with local variables? If we can do it in the former, then
what happens to the validation tests for sizes that run at the
end of the blocks --- those will also have to change to match with
tests to go along? I don't think we can do it in transformed parameters,
because that could change the number of parameters on a per-iteration
basis. So then all the discussion of matrix assignment in the doc
needs to be qualified.

Maybe we can come up with an intermediate solution, such as having
a matrix declared without any size restrictions, as in

matrix A;

And then only use this as a local variable. Or maybe in transformed
data.

Then when using that on the LHS of an array, we'd not check sizes.
That'd require a lot more code and it'd change explicitly the
sort of implicit change that Ben's arguing for. That is, we'd now
have a contrast between sized and unsized types, allowing assignment
and presumably reassignment to unsized types. That's even more code
and doc, but I think it's a much more coherent solution than modifying
the behavior of size-0 matrices.

• Bob

On Mar 26, 2015, at 8:59 AM, bgoodri notifications@github.com wrote:

If we implement true ragged arrays, that would be better because it would
work in data and parameters as well. But we've been talking about doing
that since 2011. And only Bob can do something like that with the
language.
Whereas it takes 4 lines of already-written code to assign to a 0x0
matrix
to overcome a restriction we imposed on the language for a conference
paper.

On Wed, Mar 25, 2015 at 3:13 PM, Daniel Lee notifications@github.com
wrote:

I think Bob's suggestion is what makes the most sense. I'm guessing
will
take a long time to figure out, so, let's talk about what we can do in
the
short term.

@bgoodri https://github.com/bgoodri, I still think allowing
matrix[0,0]
to mean something other than a zero-row, zero-column matrix is a bit
dangerous and could lead to model building errors that aren't clear. I
also
think having a naked matrix call is also just asking for trouble, but

it's better in my book.

I know, I didn't come up with a decent counter-example, but perhaps
something like:

data {
int N;
int M;
matrix[N,N] X1;
matrix[M,M] X2;
}
transformed data {
matrix[M,M] L_X1;
matrix[M,M] L_X2;
L_X1 <- cholesky_decompose(X1);
L_X2 <- cholesky_decompose(X2);
}
...

That'll work if M is 0. Seems odd and error-prone to me.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

[bgoodri]

If there were some way to turn checks off instead of taking them out of the
language completely, that would be okay but we haven't implemented that in
years either. To take it somewhat off track, for huge K, something like

transformed parameters {
corr_matrix[K,K] Sigma;
Sigma <- multiply_lower_tri_self_transpose(L);
}

is infeasible and may very well fail the positive definiteness check for
numerical reasons. So, the "best" way to work around this in Stan is to
duplicate code in the generated quantities block

model {
corr_matrix[K,K] Sigma;
Sigma <- multiply_lower_tri_self_transpose(L);
...
}
generated quantities {
corr_matrix[K,K] Sigma;
Sigma <- multiply_lower_tri_self_transpose(L);
}

And if the "best" thing to do is to do something that is inconsistent with
one of the central tenets of programming, then the design has a problem.
This problem could be solved by disabling the positive definiteness check
and the matrix resizing issue could be solved by disabling the size check.
But if checks like these can't be disabled, then we are punishing people
who write correct code and who want to take advantage of unimplemented
features in order to catch mistakes of people who write incorrect code.

If someone defines a corr_matrix or cov_matrix in the transformed
parameters block that is not positive definite by construction, then Stan
is drawing from an invalid posterior distribution even if the lack of
positive definiteness is caught because the priors on the parameters do not
integrate to 1 when the positive definiteness check enforces a region of
zero prior density. And there is nothing that we can do about the fact that
it is an unsound posterior and yet we do validity checks on the objects
that slow sound code down. It's fine if you want the checks to be on in a
test run to make sure that what you intended to be a positive definite
transformation is, in fact, positive definite. But it is bad to force such
checks to be on always.

Ben

On Wed, Mar 25, 2015 at 8:21 PM, Bob Carpenter notifications@github.com
wrote:

I agree about allowing declared ragged arrays anywhere. And I
fully admit I don't have this on my to-do list for any time soon.

I strongly disagree about turning off size checks at the ends of
blocks. And in general about favoring functionality over error
checking.

Maybe we can discuss at the next Stan meeting?

On Mar 26, 2015, at 10:27 AM, bgoodri notifications@github.com wrote:

If we are going to support genuine ragged arrays like

matrix[rows,cols] X[elements];

where rows and cols are vectors, then I think we should allow that
anywhere.

If it is just the original hack where a 0x0 matrix or 0x1 vector can get
reassigned to any genuine size, then it definitely can't be data or
parameters and if it can only be model or generated quantities, then that
would still be useful.

More generally, I am not a fan of limiting language features that lots of
people would find useful in order to catch bugs in user code. If we can
catch bugs like that essentially costlessly, then okay. But, such bugs
will
probably manifest themselves elsewhere and if not, too bad for that
particular user, as opposed to bad for all users. So, in this case, are
you
saying that

transformed parameters {
cov_matrix[K,K] X;
X <- crossprod(Y);
}

currently checks at the end of the block that X is, in fact, KxK. If so,
I
agree it would be hard to check if X were declared 0x0 and resized, but I
would favor dropping that check. But we could still check that it is
positive definite without regard to the size, right?

On Wed, Mar 25, 2015 at 7:05 PM, Bob Carpenter <notifications@github.com

wrote:

This project has to be a compromise. I don't like this
proposal because it messes up the clean conception of sizing
that we have everywhere. But I agree that we've been wanting
to do this forever. If I could ever stop treading water, I
could get lots of these features in. Maybe I'll try not responding
to the mailing list again and see if anyone picks up the slack.
Daniel's been volunteering.

Don't forget doc and testing. It's may be 4 lines of
code to change the assign function, but then there's testing
and mods to the doc everywhere we talk about assignment.

By the way, is this only for matrices, or also for vectors
and row vectors?

I think it'll actually require a fair amount of doc and testing.
Can we do this reassignment in the transformed data block or just
with local variables? If we can do it in the former, then
what happens to the validation tests for sizes that run at the
end of the blocks --- those will also have to change to match with
tests to go along? I don't think we can do it in transformed
parameters,
because that could change the number of parameters on a per-iteration
basis. So then all the discussion of matrix assignment in the doc
needs to be qualified.

Maybe we can come up with an intermediate solution, such as having
a matrix declared without any size restrictions, as in

matrix A;

And then only use this as a local variable. Or maybe in transformed
data.

Then when using that on the LHS of an array, we'd not check sizes.
That'd require a lot more code and it'd change explicitly the
sort of implicit change that Ben's arguing for. That is, we'd now
have a contrast between sized and unsized types, allowing assignment
and presumably reassignment to unsized types. That's even more code
and doc, but I think it's a much more coherent solution than modifying
the behavior of size-0 matrices.

• Bob

On Mar 26, 2015, at 8:59 AM, bgoodri notifications@github.com
wrote:

If we implement true ragged arrays, that would be better because it
would
work in data and parameters as well. But we've been talking about
doing
that since 2011. And only Bob can do something like that with the
language.
Whereas it takes 4 lines of already-written code to assign to a 0x0
matrix
to overcome a restriction we imposed on the language for a conference
paper.

On Wed, Mar 25, 2015 at 3:13 PM, Daniel Lee <
notifications@github.com>
wrote:

I think Bob's suggestion is what makes the most sense. I'm guessing
will
take a long time to figure out, so, let's talk about what we can
do in
the
short term.

@bgoodri https://github.com/bgoodri, I still think allowing
matrix[0,0]
to mean something other than a zero-row, zero-column matrix is a
bit
dangerous and could lead to model building errors that aren't
clear. I
also
think having a naked matrix call is also just asking for trouble,
but

it's better in my book.

I know, I didn't come up with a decent counter-example, but perhaps
something like:

data {
int N;
int M;
matrix[N,N] X1;
matrix[M,M] X2;
}
transformed data {
matrix[M,M] L_X1;
matrix[M,M] L_X2;
L_X1 <- cholesky_decompose(X1);
L_X2 <- cholesky_decompose(X2);
}
...

That'll work if M is 0. Seems odd and error-prone to me.

—
Reply to this email directly or view it on GitHub
<
https://github.com/stan-dev/stan/issues/1369#issuecomment-86178384>.

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bob-carpenter]

On Mar 26, 2015, at 12:57 PM, bgoodri notifications@github.com wrote:

If there were some way to turn checks off instead of taking them out of the
language completely, that would be okay but we haven't implemented that in
years either.

I'm OK with being able to disable all checks if someone wants to do it.
Part of the refactoring Daniel and Rob undertook was supposed to address
this or at least make it easier.

To take it somewhat off track, for huge K, something like

transformed parameters {
corr_matrix[K,K] Sigma;
Sigma <- multiply_lower_tri_self_transpose(L);
}

is infeasible and may very well fail the positive definiteness check for
numerical reasons. So, the "best" way to work around this in Stan is to
duplicate code in the generated quantities block

model {
corr_matrix[K,K] Sigma;
Sigma <- multiply_lower_tri_self_transpose(L);
...
}

You can't declare a corr_matrix in the model, but you could
do this with matrix[K]. Local variables don't get checked
for constraints, which is why you can't define constraints on them.

So if you define Sigma as a matrix[K,K] in the transformed parameters
block, it won't do the checks either.

generated quantities {
corr_matrix[K,K] Sigma;
Sigma <- multiply_lower_tri_self_transpose(L);
}

And if the "best" thing to do is to do something that is inconsistent with
one of the central tenets of programming, then the design has a problem.

Agreed. This seems like a problem with our enforcing positive
definiteness.

This problem could be solved by disabling the positive definiteness check
and the matrix resizing issue could be solved by disabling the size check.

Where would it be OK to use a correlation matrix that doesn't pass
positive-definiteness checks? Shouldn't it then just be declared
as a matrix?

But if checks like these can't be disabled, then we are punishing people
who write correct code and who want to take advantage of unimplemented
features in order to catch mistakes of people who write incorrect code.

What I'm missing is how the code is correct if they're
producing non-positive definite matrices and assigning them to
correlation matrix variables.

Is the problem just that we need a different test?

If someone defines a corr_matrix or cov_matrix in the transformed
parameters block that is not positive definite by construction, then Stan
is drawing from an invalid posterior distribution even if the lack of
positive definiteness is caught because the priors on the parameters do not
integrate to 1 when the positive definiteness check enforces a region of
zero prior density. And there is nothing that we can do about the fact that
it is an unsound posterior and yet we do validity checks on the objects
that slow sound code down. It's fine if you want the checks to be on in a
test run to make sure that what you intended to be a positive definite
transformation is, in fact, positive definite. But it is bad to force such
checks to be on always.

I'm confused, because this seems to be arguing that we should enforce
the checks to prevent drawing from an invalid posterior.

We could go further and just terminate the run if we get a "divergent"
iteration in Michael sense and thus force users to step down the step
sizes.

• Bob

[syclik]

Ben, if you want us to turn off all checks, that isn't the hardest thing to
plumb all the way through. The issue is that I think you underestimate how
much Eigen and Boost likes to swallow errors. For example, the llt
transform will take nans, infs, pretty much anything, and give back an "llt
transform" that has no indicators that it's broken. I don't think they doc
it as a precondition, but the transforms are obviously taking in bad
matrices and then pushing out things that everything else thinks is good.

If you're ok with that, we can go ahead and find a way to disable all
checks. It just makes me really, really nervous when we know that
something, at runtime, will mask something bad. If you really want it in,
create a new issue, give us some hints as to what you want. All checks off?
A subset of checks off? matrix checks off? bounds? hopefully you get the
point.

What I'm missing is how the code is correct if they're
producing non-positive definite matrices and assigning them to
correlation matrix variables.

Is the problem just that we need a different test?

I have the same question: what should we do?

[bgoodri]

On Wed, Mar 25, 2015 at 10:21 PM, Bob Carpenter notifications@github.com
wrote:

And if the "best" thing to do is to do something that is inconsistent
with
one of the central tenets of programming, then the design has a problem.

Agreed. This seems like a problem with our enforcing positive
definiteness.

This problem could be solved by disabling the positive definiteness check
and the matrix resizing issue could be solved by disabling the size
check.

Where would it be OK to use a correlation matrix that doesn't pass
positive-definiteness checks? Shouldn't it then just be declared
as a matrix?

It is a dilemma that is self-inflicted by the user who writes
transformations that are not consistent with the constraints by
construction.

• If someone passes an indefinite symmetric matrix to a function like
  multi_normal, then any number multi_normal() outputs is not sound
  (hopefully it would throw an exception or return NaN).
• If someone passes only the symmetric matrices that pass the positive
  definiteness check to a function like multi_normal, then multi_normal won't
  throw an exception or return NaN but the draws are not sound. Basically,
  the user should be truncating the relevant priors by dividing them by the
  non-constant 1 - Pr(indefinite) but there is no way to do that truncation
  in Stan.

So, the only way to get sound posterior draws is to only use
transformations that are positive definite with probability one. And if do
that, then you don't want to check for positive definiteness at runtime in
production code both because it is slower and also because you can get
positive definite matrices that are numerically indefinite.

For positive definiteness, there are ways of circumventing the check by
just declaring it as a plain square matrix. But the larger point I am
trying to make is that requiring everything to be checked for validity
hurts people who write valid code in order to prevent people who write
invalid code from hurting themselves. In the original case of size checking
matrices, there is virtually no runtime cost but the cost is that we can't
resize them and we can't implement C++ or user-defined functions that
return multiple matrices of different sizes. So Sebastian can't do things
like

transformed parameters {
matrix my_stuff[P];
my_stuff <- make_stuff(lots, of, arguments); // for a user-defined
make_stuff
}

and export that single make_stuff function to R to use in posterior
predictive simulations. But it is really the principle that everyone has to
suffer in order to catch someone else's typos that bothers me.

Ben

[syclik]

@bgoodri, do you have an example of a positive definite matrix that is numerically indefinite? (I'm trying to search online for one, but can't find one easily)

[bgoodri]

You can make one in Stan with something like

L <- lkj_cholesky_corr_rng(K);
d <- append_col(rep_vector(1, K - 1), 1e-20);
Sigma <- multiply_lower_tri_self_transpose(diag_pre_multiply(d, L));

On Wed, Mar 25, 2015 at 11:01 PM, Daniel Lee notifications@github.com
wrote:

@bgoodri https://github.com/bgoodri, do you have an example of a
positive definite matrix that is numerically indefinite? (I'm trying to
search online for one, but can't find one easily)

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[syclik]

I've managed to get an example running and I see it's failing, but is that an artifact of going to disk and back? Could we have issues with: - multiply_lower_tri_self_transpose? - diag_pre_multiply? or do you really believe those are solid enough to not produce numerically unstable results? There are too many moving parts, so I was hoping you could provide an example that I could look at. (and hopefully it's smaller than the one I generated?)

[bgoodri]

Well, it is that multiply_lower_tri_self_transpose produces numerically
unstable results. Even though it maximizes the precision, it can't preserve
enough precision for a Cholesky factorization to come out with positive
diagonals. Just

cov_matrix[2] Sigma;
Sigma[1,1] <- 1;
Sigma[1,2] <- 1 - 1e-17;
Sigma[2,1] <- Sigma[1,2]
Sigma[2,2] <- 1;

will underflow the off-diagonal element and result in a numerically
singular matrix.

On Thu, Mar 26, 2015 at 12:45 AM, Daniel Lee notifications@github.com
wrote:

I've managed to get an example running and I see it's failing, but is that
an artifact of going to disk and back?

Could we have issues with:

• multiply_lower_tri_self_transpose?
• diag_pre_multiply?

or do you really believe those are solid enough to not produce numerically
unstable results?

There are too many moving parts, so I was hoping you could provide an
example that I could look at. (and hopefully it's smaller than the one I
generated?)

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bgoodri]

How about we add this Eigen::Matrix constructor by adding this to MatrixAddons.h

EIGEN_STRONG_INLINE Matrix(const std::vector<int> &x)
  : Matrix(x[0],x[1]) // only works for C++11 but we could figure something out
{
}

Then amend the Stan language to accept

int<lower=1> indices[J,2];
...
matrix[indices] X[J];

and that parses using the range constructor
http://www.cplusplus.com/reference/vector/vector/vector/
to

std::vector<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > X(indices.begin(), indices.end());

[syclik]

I think that looks great. Think we'd be able to declare matrix sizes using
that format? So something like:

int indices[2];
matrix[indices] Y;

Or would this just be available an array of matrices so the equivalent
would be:
int indices[1,2];
matrix[indices] Y[1];

On Fri, Mar 27, 2015 at 5:51 PM, bgoodri notifications@github.com wrote:

How about we add this Eigen::Matrix constructor by adding this to
MatrixAddons.h

EIGEN_STRONG_INLINE Matrix(const std::vector &x)
: Matrix(x[0],x[1]) // only works for C++11 but we could figure something out
{
}

Then amend the Stan language to accept

int<lower=1> indices[J,2];
...
matrix[indices] X[J];

and that parses using the range constructor
http://www.cplusplus.com/reference/vector/vector/vector/
to

std::vectorEigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic > X(indices.begin(), indices.end());

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bgoodri]

I guess we could also support
matrix[indices] Y; // no array

On Fri, Mar 27, 2015 at 6:05 PM, Daniel Lee notifications@github.com
wrote:

I think that looks great. Think we'd be able to declare matrix sizes using
that format? So something like:

int indices[2];
matrix[indices] Y;

Or would this just be available an array of matrices so the equivalent
would be:
int indices[1,2];
matrix[indices] Y[1];

On Fri, Mar 27, 2015 at 5:51 PM, bgoodri notifications@github.com wrote:

How about we add this Eigen::Matrix constructor by adding this to
MatrixAddons.h

EIGEN_STRONG_INLINE Matrix(const std::vector &x)
: Matrix(x[0],x[1]) // only works for C++11 but we could figure
something out
{
}

Then amend the Stan language to accept

int<lower=1> indices[J,2];
...
matrix[indices] X[J];

and that parses using the range constructor
http://www.cplusplus.com/reference/vector/vector/vector/
to

std::vectorEigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic >
X(indices.begin(), indices.end());

—
Reply to this email directly or view it on GitHub
#1369 (comment).

—
Reply to this email directly or view it on GitHub
#1369 (comment).

[bob-carpenter]

I like the basic idea --- it's along the lines of what I've been
proposing for more general declarations. We could do similar things
with constraints so that not every element of a vector gets the same
bounds.

I don't quite understand the details, though.

I think I'd prefer to not use 2D arrays as implicit codings.
How about something like:

int J;
int row_sizes[J];
int col_sizes[J];
matrix[row_sizes,col_sizes] x[J];

I'll have to think through whatever else is going to come up from this.

I think we should just code this up directly in Stan rather than trying to
use MatrixAddons.h --- we finally got rid of it in favor of templates
for scalar and index types. The problem is that it induces include
order complexity that nobody could keep straight when writing matrix
functions.

• Bob

P.S. I'm struggling to just tread water with keeping existing things
working and follow all the requests. I have no idea how I'm going to
free up time to write new code.

On Mar 28, 2015, at 8:51 AM, bgoodri notifications@github.com wrote:

How about we add this Eigen::Matrix constructor by adding this to MatrixAddons.h

EIGEN_STRONG_INLINE Matrix(const std::vector &x)
: Matrix(x[0],x[1]) // only works for C++11 but we could figure something out
{
}

Then amend the Stan language to accept

int<lower=1> indices[J,2];
...
matrix[indices] X[J];

and that parses using the range constructor
http://www.cplusplus.com/reference/vector/vector/vector/
to

std::vectorEigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic > X(indices.begin(), indices.end());

—
Reply to this email directly or view it on GitHub.

[bgoodri]

On Fri, Mar 27, 2015 at 7:31 PM, Bob Carpenter notifications@github.com
wrote:

How about something like:

int J;
int row_sizes[J];
int col_sizes[J];
matrix[row_sizes,col_sizes] x[J];

Fine with me, but ...

P.S. I'm struggling to just tread water with keeping existing things
working and follow all the requests. I have no idea how I'm going to
free up time to write new code.

No one else can make non-trivial changes to the Stan language, so I am only
proposing ideas that require small changes, like changing assign.hpp or
iterating over a vector of two indices. I don't know how to suggest parsing

int J;
int row_sizes[J];
int col_sizes[J];
matrix[row_sizes,col_sizes] x[J];

Ben

[bob-carpenter]

On Mar 28, 2015, at 10:38 AM, bgoodri notifications@github.com wrote:

On Fri, Mar 27, 2015 at 7:31 PM, Bob Carpenter notifications@github.com
wrote:

How about something like:

int J;
int row_sizes[J];
int col_sizes[J];
matrix[row_sizes,col_sizes] x[J];

Fine with me, but ...

P.S. I'm struggling to just tread water with keeping existing things
working and follow all the requests. I have no idea how I'm going to
free up time to write new code.

No one else can make non-trivial changes to the Stan language,

Daniel says he understands it enough to do so, but I don't know if
he's made major changes. Marcus added the unit_length data type.
Mitzi added the exceptions and the operator^.

I think we really need to agree on priorities. I just feel like
I'm being pulled in like ten directions at once and I can't scale.

so I am only
proposing ideas that require small changes, like changing assign.hpp or
iterating over a vector of two indices. I don't know how to suggest parsing

int J;
int row_sizes[J];
int col_sizes[J];
matrix[row_sizes,col_sizes] x[J];

But your suggested change would also require changing the parser and
code generator and all the validating code.

• Bob

[bgoodri]

On Fri, Mar 27, 2015 at 7:49 PM, Bob Carpenter notifications@github.com
wrote:

But your suggested change would also require changing the parser and
code generator and all the validating code.

I don't know how many files would have to change, but I assumed the
validating code for an array is already something like

for (j in 1:J) validate(X[j]);

in which case it seems like it should work the same whether the elements of
the array have homogenous or heterogenous sizes. I suppose that if
validate() also checks whether X[j] has the size that X[j] was declared
with, then that would have to be amended, but I don't even know how you
could construct a matrix currently that would fail a size check.

Ragged arrays have been one of the most requested features since 2011. They
are necessary to do things like

• Matrix decompositions that return factors of different sizes like QR
  without calling qr_Q() and qr_R()
• Kronecker products of an arbitrary number of arbitrarily sized
  covariance matrices that are needed by Seth / Aki / Peter
• Writing all of transformed parameters as one call to a user-defined
  function that returns a bunch of matrices that can be exported to R /
  tested / called after the sampling is done
• Doing multilevel models without a crazy amount of copying or index
  fiddling

and probably a bunch of other use cases. If anyone else could have
implemented them, they probably would have.

Ben

[bob-carpenter]

On Mar 28, 2015, at 11:00 AM, bgoodri notifications@github.com wrote:

On Fri, Mar 27, 2015 at 7:49 PM, Bob Carpenter notifications@github.com
wrote:

But your suggested change would also require changing the parser and
code generator and all the validating code.

I don't know how many files would have to change, but I assumed the
validating code for an array is already something like

for (j in 1:J) validate(X[j]);

in which case it seems like it should work the same whether the elements of
the array have homogenous or heterogenous sizes.

I think you're right here, but I'd have to look and we'd need
tests around all of it to make sure.

I suppose that if
validate() also checks whether X[j] has the size that X[j] was declared
with, then that would have to be amended, but I don't even know how you
could construct a matrix currently that would fail a size check.

Right --- we don't check that because we don't allow resizing
assignments. So the assignment checking would have to change, but
I believe that's what you're proposing.

Ragged arrays have been one of the most requested features since 2011.

Stan 1.0 was only released August 2012!

Adding line numbers to error messages was also on the list
since before Stan 1.0, and I just got to it. I feel pretty strongly
that we should prioritize making things that work solid rather than
adding new features.

I did add user-defined functions last year and all the language
stuff around ODEs, and those have also been on the list from the
get go. They then brought up a ton of issues which needed to be
resolved.

If I could make more time to do things, I would. I really feel
bad I can't get to everyone's (or most of the time anyone's) feature
requests.

Right now, my priorities are two grant proposals and finishing the
autodiff paper. I have completely ignored the feature requests
from Michael and Andrew for extensions to the language to allow
stochastic VB, whatever we need for MML, and adiabatic sampling.
We really need to get higher-order autodiff done for MLE.

See my recently added longer-term to-do list on the stan-dev Wiki.
Let's talk about priorities at the next Stan meeting.

I'm going to start by not answering any more issues on stan-users
and see if anyone else steps up to the plate. I can just send
an "out-of-office" auto-reply that says we're too busy to help if
they don't get answered in a day or two. I think this'll seriously
hurt our reputation as having great support, but maybe new features
are more important.

I also spend a lot of time trying to stay on top of issues like
command interfaces and API refactoring and lint checking, mainly
because I disagree with some of the proposals I see. Maybe I just
need to learn to let go. I obviously can't scale.

Having said that, I really don't want to keep throwing down hacks.
It makes things impossible to maintain going forward both in terms
of docs and testing and future coding.

They
are necessary to do things like

• Matrix decompositions that return factors of different sizes like QR
  without calling qr_Q() and qr_R()
• Kronecker products of an arbitrary number of arbitrarily sized
  covariance matrices that are needed by Seth / Aki / Peter
• Writing all of transformed parameters as one call to a user-defined
  function that returns a bunch of matrices that can be exported to R /
  tested / called after the sampling is done
• Doing multilevel models without a crazy amount of copying or index
  fiddling

and probably a bunch of other use cases. If anyone else could have
implemented them, they probably would have.

Mitzi's saying she'll try to take it on. First, we need a design
for how it'll work.

• Bob

[bgoodri]

On Fri, Mar 27, 2015 at 8:25 PM, Bob Carpenter notifications@github.com
wrote:

Ragged arrays have been one of the most requested features since 2011.

Stan 1.0 was only released August 2012!

It was on the TODO list well before 1.0

https://groups.google.com/d/msg/stan-dev/TrvVDfQ9uMc/Yt9n0dbeSHwJ

If Mitzi can get the ball rolling, that could make the difference.

Ben