# Scalar, vector, and matrix views

#1

Is there some reason we can’t replace our basic views with
something like the following. I know we had some writable vector
views, but can those be put into a different class?

• Bob

// wrap (scalar OR std::vector OR Eigen vector OR Eigen row
// vector) to present sequence of scalars
template
struct ScalSeqView { };

template
struct ScalSeqView {
const T& x_;
explicit SeqView(const T& x) : x_(x) { }
const T& operator[](int n) const { return x_; }
};

template
struct ScalSeqView<std::vector > {
const vector& xs_;
explicit SeqView(const vector& xs) : xs_(xs) { }
const T& operator[](int n) const { return xs_[n]; }
};

template
struct ScalSeqView<Matrix<T, 1, -1> > {
const Matrix<T, 1, -1>& xs_;
explicit SeqView(const Matrix<T, 1, -1>& xs) : xs_(xs) { }
const T& operator[](int n) const { return xs_[n]; }
};

template
struct ScalSeqView<Matrix<T, -1, 1> > {
const Matrix<T, -1, 1>& xs_;
explicit SeqView(const Matrix<T, -1, 1>& xs) : xs_(xs) { }
const T& operator[](int n) const { return xs_[n]; }
};

// wrap Eigen vector or std::vector of Eigen vector
// to produce Eigen vector; also works for Eigen row vectors
// and Eigen matrices
template <typename S, R, C>
struct MatrixSeqView { };

template <typename T, R, C>
struct MatrixSeqView<Matrix<T, R, C> > {
const Matrix<T, R, C>& x_;
explicit MatrixSeqView(const Matrix<T, R, C>& x) { x_(x); }
const Matrix<T, R, C>& operator[](int n) const { return x_; }
};

template <typename T, R, C>
MatrixSeqView<vector<Matrix<T, R, C> >, R, C>> {
const vector<Matrix<T, R, C> >& xs_;
MatrixSeqView(const vector<Matrix<T, R, C> >& xs) { xs_(xs); }
const Matrix<T, -1, -1>& operator[](int n) const { return xs_[n]; }
};

#2

This seems like it drops the “broadcasting” behavior we talked about, where no matter what index is passed in the original `T` was returned, right?

Given a matrix, is it supposed to iterate through it scalar by scalar?

#3

This was a very difficult place to have you start looking
at Stan as it’s poorly written and named code that’s tied up
in some of our most complicated template metaprogramming. Maybe we
can back off and give you something simpler to work on first?
But let me answer.

There are three different wrappers we need:

(a) return scalars [VectorView]
(b) return Eigen vectors
© return Eigen matrices [VectorViewMvt]

I didn’t write most of this code, so I’m not exactly sure how
everything’s plumbed in. But we need the above eventually to
"vectorize" scalars, vectors, and matrices.

I’m not sure what you mean by “given a matrix”. Which of
(a–c)? I’m not sure we ever use VectorView of Eigen matrices
right now, though I think it’s implemented for that. I’m
not sure how we’re doing (b) now, but it may just be ad-hoc
code in some of the distributions with vector parameters (like
multivariate normal, which also has a matrix parameter).

• Bob

#4

Gotcha. I am happy to work on whatever you guys find useful and you think would make a good starting point. As a straw man proposal, how about I add tests for the currently used and correct behavior in `VectorViewMvt` and `scalar_type_pre` so we can close #93 with it, and leave the refactoring and bug in separate issues?

And then do you have any suggestions for next issues? I took a look and mostly have no idea what I’m looking at, but a couple ideas stood out as maybe easy enough:

#5

That’d be great. One thing we want to do is expose
all of our transforms. The problem we’re having with
a lot of these issues is that we’ve already done the
simple ones. A great project would be to expose the
transforms, but that has tentacles all over the code and
we want to combine it with refactoring them and vectorizing
them. And I think Thel may already be working on it, but
I’m not sure. Thel, are you listening?

• Bob

#6

Would be a good separate thread to come up with projects that are do-able and priorities without diving into all the innards of Stan’s template programming.

#7

@sakrejda, please don’t hijack threads. And to answer your question / non-question: yes. Please feel free to start a new thread.

#8

@thel, there was a question directed at you.

#9

I don’t know who originally wrote it. I’m guessing it’s missing doc. You’re down in the guts of Stan where we haven’t actually touched a lot of it in a long time. We really need to go back and doc everything and test it. Clearly, it’s hard to move forward without knowing what the current stuff does.

I’m in favor of renaming things as Bob wrote in the original post.

#10

I’m watching the thread and just wanted to chime in that nothing I’m doing is likely to obstruct work on this stuff. I’d be happy to get into rewriting some of it as it goes, please tag me in or add me as a watcher on github.

#11

Rob Trangucci wrote off list:

There is a class we added about 2 years ago called container_view that mirrors the `VectorView`. It was intended to be used in a multi-variate version of `OperandsAndPartials`. A multi-variate version of `OperandsAndPartials` will be very useful for writing a version of `multi_normal_cholesky` that has analytical gradients.

I’d like to leave the `container_view` in the math library, but it should probably be refactored to mirror the changes you’ll make to the `VectorView` class.

#12

I didn’t see `container_view` used anywhere, just defined. Is it used somewhere outside of `stan-dev/math`?

`VectorViewMvt` is the multivariate version of `VectorView`. We’ll be consolidating that and `VectorView` into these three classes (and maybe a fourth for `row_vector_seq_view` if necessary):

• `scalar_seq_view`: `real, real[], vector, row_vector`
`const real& operator[](int)`

• `vector_seq_view`: `vector, vector[]`
`const vector& operator[](int)`

• `matrix_seq_view`: `matrix, matrix[]`
`const matrix& operator[](int)`

What Rob’s talking about is here:

``````stan/math/prim/scal/meta/container_view.hpp
stan/math/prim/arr/meta/container_view.hpp
stan/math/prim/mat/meta/container_view.hpp
``````

I think I may need some help understanding the motivation of Rob’s `container_view`. What puzzles me is that I don’t see it used anywhere on the `develop` branch. Did this get put in on spec for some future work?

If we need it, I think we can improve the design to use specializations for the various container types with reference member variables; that way we won’t even need a `Map`. I don’t see what the matrix version of this gives you over just using a `Map` directly, though.

Also, is it intended to be mutable? It will be the way it’s designed, since that `Map` (in the matrix version) will hold a reference to the original memory.

I want to define a separate set of mutable sequence views. We can then replace the builders used in `OperandsAndPartials` with those.

#13

The intention was for `container_view` to offer a matrix/array of matrices generalization one of the 3 behaviors in VectorView: allowing a mutable scalar view of an array of doubles. `container_view` provides a way to return a Vector, Matrix, or a `std::vector` of vectors/matrices views of arrays of doubles. This would allow us to build a version of `operands_and_partials` where the gradients could have matrix structure, intended for use something like `multi_normal_cholesky`. We haven’t built/modified multivariate `operands_and_partials` yet, but there is a clear need for something like this to simplify the coding of multivariate densities.

`VectorViewMvt` only does array broadcasting for vectors and matrices because it wasn’t intended to be integrated into a multivariate version of `operands_and_partials`. Perhaps `container_view` isn’t the right name for the class, and should be changed to whatever you call the mutable array version of `VectorView`.

#14

That’s a different application than `T_seq_view`, the point of which is to allow a `T` or a container of `T` to be viewed as a sequence of `T`.

What I don’t understand is why you need another class. Can’t you just use the `Map` directly? Or are you thinking `OperandsAndPartials` will hold pointers to these map holders which will point into their raw memory?

#15

That’s fine; we can rename `container_view` whatever you’ll be renaming the mutable version of VectorView on lines 11 - 24: https://github.com/stan-dev/math/blob/1f2794953e5199717afd87bf40dbbe1f7165ce20/stan/math/prim/mat/meta/VectorView.hpp

We could do this. It just isn’t how `OperandsAndPartials` is currently structured. Will you be rewriting https://github.com/stan-dev/math/blob/develop/stan/math/rev/scal/meta/OperandsAndPartials.hpp to not rely on VectorView? If so, then we can generalize `OperandsAndPartials` for multivariate operands and bring the `container_view` functionality into `OperandsAndPartials`.

#16

Yes, completely getting rid of `vector_view`. The `OperandsAndPartials` class has two unusual use cases for it, so we’ll do those last. My goal is to tease apart the usages into separate classes rather than have these things brimming with template parameters and indirection in attempt to perform multiple functions. So that’s not going to be in the first pull request, I don’t think.

#17

Great! So maybe it would make sense to try to generalize `OperandsAndPartials` after vector_view. We can rename container_view to something more informative for use in the general version. The general `OperandsAndPartials` will really simplify writing gradients with respect to matrix-variate parameters.

#18

I’m taking a look at `OperandsAndPartials` generally (as part of looking at how it uses `VectorView`) and I have two questions:

1. Are we okay with removing the guards around the `throw_if_accessed` behavior in `VectorView`? This would mean when we end up using `OperandsAndPartials` for constants we don’t have the protection against the programmer accidentally accessing its elements.
2. Above we talk a bit about generalizing `OperandsAndPartials` to the multivariate case. Looking at the implementation of `normal_lpdf` vs `multi_normal_lpdf` they use a similar number of lines and I don’t think I have enough underlying math to totally understand how a multivariate version would improve the code in the multinormal case. @rtrangucci would you want to sit down and go over what you’re thinking for `OperandsAndPartials` and maybe explain a bit how it’ll be used? I think that’d help me a bunch in this area, though I could also just leave this part of it for later or someone else.

Thanks,
Sean

#19

Yes, absolutely. Bob convinced me a while ago that we don’t have to really defend against all programmers. It was originally set up that way because it’s easy to instantiate things that run over memory without realizing it. There were a couple bugs that we caught with that sort of behavior, but I think we’re good now. Especially if we’re able to doc it properly.

For the simple case, we would be looking at vectorizing all the multivariate distributions by using a view, not necessarily with operands and partials. We would let autodiff work its magic. The next level up would be to build up the operands and partials using `double` matrices, but that’s a lot more effort.

#20

We have this currently with the VectorViewMvt. I think we should do the multivariate version of `OperandsAndPartials`, and I do have time to meet with @seantalts this week. Let’s meet on Thursday after the meeting.