Deducing Return type for new matrix when no inputs are matrix<var> or var<matrix>

[SteveBronder]

Description

For the var<Matrix> our current rules are that if any of the inputs are var<Matrix> the output will also be var<Matrix>. But there's several functions like pow() which can take in non-var matrix types and produce a var<Matrix>. i.e. like the signature

// elementwise pow
Matrix<var> pow(Matrix<double>, var);
Matrix<var> gp_exp_quad_cov(std::vector<double>, var, var)

For these I think there are two schemes.

1. A template parameter Ret like with rep_vector() for the var specialization with a default of an Eigen::Matrix. Then in the compiler we can change Ret to a var_value<Eigen::Matrix> by passing a template value.

This gets weird because idk how to define Ret when a function can take in a matrix, vector, or row_vector.

2. A template parameter is_var_matrix with a default value of false.

Little weirder but a lot easier. Pow for example would be something like the below, where conditional_var_matrix_t takes in a bool and a type and if the bool is true returns a var<Matrix>.

template <bool is_var_matrix = false, typename Mat, require_matrix_t<Mat>* = nullptr>
auto pow(const Mat& base, var exponent) {
 // at some point in the function have a type trait to deduce whether to return a var<Matrix> or not
 conditional_var_matrix_t<is_var_matrix, Mat> ret = ...
}

and for gp_exp_quad_cov I think we actually need to functions, so in that case we would do SFINAE off of is_var_matrix

template <bool is_var_matrix = false, typename T_x, typename T_sigma, require_st_arithmetic<T_x>* = nullptr,
          require_stan_scalar_t<T_sigma>* = nullptr,
         // Checking for return type 
         require_t<bool_constant<!is_var_matrix>>* = nullptr>
inline Eigen::Matrix<var, -1, -1> gp_exp_quad_cov(const std::vector<T_x>& x,
                                                  const T_sigma sigma,
                                                  const var length_scale) {
  // at some point in the function 
  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
}

template <bool is_var_matrix = false, typename T_x, typename T_sigma, require_st_arithmetic<T_x>* = nullptr,
          require_stan_scalar_t<T_sigma>* = nullptr,
         // Checking for return type 
         require_t<bool_constant<is_var_matrix>>* = nullptr>
inline Eigen::Matrix<var, -1, -1> gp_exp_quad_cov(const std::vector<T_x>& x,
                                                  const T_sigma sigma,
                                                  const var length_scale) {
  // at some point in the function 
  var_value<Eigen::Matrix<double, -1, -1>> cov(Eigen::Matrix<double, -1, -1>(x_size, x_size));
}

Current Version:

v4.1.0

[t4c1]

This gets weird because idk how to define Ret when a function can take in a matrix, vector, or row_vector.

You would have to make separate overloads for these three options.

I prefer the option 1, because the same pattern can be used for deciding whether to use OpenCL overloads.

[SteveBronder]

Yes as I'm working over the ocaml in the compiler I think 1 can work nicely

[SteveBronder]

The list we need to target is essentially

1. Anything that takes in a scalar and outputs a matrix
2. Anything that takes in a combination of scalars and matrices (for f(var scalar, data matrix))