Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 37 additions & 19 deletions src/stan/model/indexing/assign.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@ namespace model {
* index_max - index from 1:max
* index_min_max - index from min:max
* nil_index_list - no-op
* Ranges are empty when the lower bound exceeds the upper bound, including
* min > N for index_min and max < 1 for index_max.
* The order of the overloads are
* vector / row_vector:
* - all index overloads
Expand Down Expand Up @@ -169,11 +171,16 @@ template <typename Vec1, typename Vec2,
require_all_vector_t<Vec1, Vec2>* = nullptr,
require_all_not_std_vector_t<Vec1, Vec2>* = nullptr>
inline void assign(Vec1&& x, const Vec2& y, const char* name, index_min idx) {
stan::math::check_range("vector[min] assign", name, x.size(), idx.min_);
stan::math::check_size_match("vector[min] assign", name,
x.size() - idx.min_ + 1, "right hand side",
y.size());
internal::assign_impl(x.tail(x.size() - idx.min_ + 1), y, name);
if (likely(idx.min_ <= x.size())) {
stan::math::check_range("vector[min] assign", name, x.size(), idx.min_);
stan::math::check_size_match("vector[min] assign", name,
x.size() - idx.min_ + 1, "right hand side",
y.size());
internal::assign_impl(x.tail(x.size() - idx.min_ + 1), y, name);
} else {
stan::math::check_size_match("vector[min > size] assign", name, 0,
"right hand side", y.size());
}
}

/**
Expand Down Expand Up @@ -329,13 +336,20 @@ template <typename Mat1, typename Mat2,
require_dense_dynamic_t<Mat1>* = nullptr,
require_matrix_t<Mat2>* = nullptr>
inline void assign(Mat1&& x, const Mat2& y, const char* name, index_min idx) {
const auto row_size = x.rows() - (idx.min_ - 1);
stan::math::check_range("matrix[min] assign row", name, x.rows(), idx.min_);
stan::math::check_size_match("matrix[min] assign rows", name, row_size,
"right hand side rows", y.rows());
stan::math::check_size_match("matrix[min] assign columns", name, x.cols(),
"right hand side columns", y.cols());
internal::assign_impl(x.bottomRows(row_size), y, name);
if (likely(idx.min_ <= x.rows())) {
stan::math::check_range("matrix[min] assign row", name, x.rows(), idx.min_);
const auto row_size = x.rows() - idx.min_ + 1;
stan::math::check_size_match("matrix[min] assign rows", name, row_size,
"right hand side rows", y.rows());
stan::math::check_size_match("matrix[min] assign columns", name, x.cols(),
"right hand side columns", y.cols());
internal::assign_impl(x.bottomRows(row_size), y, name);
} else {
stan::math::check_size_match("matrix[min > rows] assign rows", name, 0,
"right hand side rows", y.rows());
stan::math::check_size_match("matrix[min] assign columns", name, x.cols(),
"right hand side columns", y.cols());
}
}

/**
Expand Down Expand Up @@ -685,13 +699,17 @@ template <typename Mat1, typename Mat2, typename Idx,
require_dense_dynamic_t<Mat1>* = nullptr>
inline void assign(Mat1&& x, const Mat2& y, const char* name,
const Idx& row_idx, index_min col_idx) {
const auto start_col = col_idx.min_ - 1;
const auto col_size = x.cols() - start_col;
stan::math::check_range("matrix[..., min] assign column", name, x.cols(),
col_idx.min_);
stan::math::check_size_match("matrix[..., min] assign columns", name,
col_size, "right hand side columns", y.cols());
assign(x.rightCols(col_size), y, name, row_idx);
if (likely(col_idx.min_ <= x.cols())) {
stan::math::check_range("matrix[..., min] assign column", name, x.cols(),
col_idx.min_);
const auto col_size = x.cols() - col_idx.min_ + 1;
stan::math::check_size_match("matrix[..., min] assign columns", name,
col_size, "right hand side columns", y.cols());
assign(x.rightCols(col_size), y, name, row_idx);
} else {
stan::math::check_size_match("matrix[..., min > cols] assign columns", name,
0, "right hand side columns", y.cols());
}
}

/**
Expand Down
31 changes: 22 additions & 9 deletions src/stan/model/indexing/rvalue.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ namespace model {
* index_max - index from 1:max
* index_min_max - index from min:max
* nil_index_list - no-op
* Ranges are empty when the lower bound exceeds the upper bound, including
* min > N for index_min and max < 1 for index_max.
* The order of the overloads are
* vector / row_vector:
* - all index overloads
Expand Down Expand Up @@ -211,8 +213,12 @@ inline auto rvalue(Vec&& v, const char* name, index_min_max idx) {
template <typename Vec, require_vector_t<Vec>* = nullptr,
require_not_std_vector_t<Vec>* = nullptr>
inline auto rvalue(Vec&& x, const char* name, index_min idx) {
stan::math::check_range("vector[min] indexing", name, x.size(), idx.min_);
return x.tail(x.size() - idx.min_ + 1);
if (idx.min_ <= x.size()) {
stan::math::check_range("vector[min] indexing", name, x.size(), idx.min_);
Comment thread
WardBrian marked this conversation as resolved.
return x.tail(x.size() - idx.min_ + 1);
} else {
return x.tail(0);
}
}

/**
Expand Down Expand Up @@ -301,9 +307,12 @@ inline auto rvalue(EigMat&& x, const char* name, MultiIndex&& idx) {
*/
template <typename Mat, require_dense_dynamic_t<Mat>* = nullptr>
inline auto rvalue(Mat&& x, const char* name, index_min idx) {
const auto row_size = x.rows() - (idx.min_ - 1);
math::check_range("matrix[min] row indexing", name, x.rows(), idx.min_);
return x.bottomRows(row_size);
if (idx.min_ <= x.rows()) {
math::check_range("matrix[min] row indexing", name, x.rows(), idx.min_);
return x.bottomRows(x.rows() - idx.min_ + 1);
} else {
return x.bottomRows(0);
}
}

/**
Expand Down Expand Up @@ -637,10 +646,14 @@ inline auto rvalue(Mat&& x, const char* name, Idx&& row_idx,
template <typename Mat, typename Idx, require_dense_dynamic_t<Mat>* = nullptr>
inline auto rvalue(Mat&& x, const char* name, Idx&& row_idx,
index_min col_idx) {
const Eigen::Index col_size = x.cols() - (col_idx.min_ - 1);
math::check_range("matrix[..., min] column indexing", name, x.cols(),
col_idx.min_);
return rvalue(x.rightCols(col_size), name, std::forward<Idx>(row_idx));
if (col_idx.min_ <= x.cols()) {
math::check_range("matrix[..., min] column indexing", name, x.cols(),
col_idx.min_);
const Eigen::Index col_size = x.cols() - col_idx.min_ + 1;
return rvalue(x.rightCols(col_size), name, std::forward<Idx>(row_idx));
} else {
return rvalue(x.rightCols(0), name, std::forward<Idx>(row_idx));
}
}

/**
Expand Down
2 changes: 1 addition & 1 deletion src/stan/model/indexing/rvalue_index_size.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ inline int rvalue_index_size(const index_omni& idx, int size) noexcept {
* @return Size of result.
*/
inline int rvalue_index_size(const index_min& idx, int size) noexcept {
return size - idx.min_ + 1;
return (idx.min_ > size) ? 0 : (size - idx.min_ + 1);
}

/**
Expand Down
4 changes: 2 additions & 2 deletions src/test/unit/model/indexing/assign_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -615,7 +615,7 @@ TEST(ModelIndexing, lvalueMatrixMultiMulti) {
EXPECT_FLOAT_EQ(y(1, 2), x(2, 3));

test_throw(x, y, index_min_max(2, 3), index_min(0));
test_throw(x, y, index_min_max(2, 3), index_min(10));
test_throw_ia(x, y, index_min_max(2, 3), index_min(10));
test_throw_ia(x, y, index_min_max(1, 3), index_min(2));

x << 0.0, 0.1, 0.2, 0.3, 1.0, 1.1, 1.2, 1.3, 2.0, 2.1, 2.2, 2.3;
Expand Down Expand Up @@ -1113,7 +1113,7 @@ TEST(model_indexing, assign_densemat_densemat_min_max_index_min_index) {
EXPECT_FLOAT_EQ(y(1, 2), x(2, 3));

test_throw(x, y, index_min_max(2, 3), index_min(0));
test_throw(x, y, index_min_max(2, 3), index_min(10));
test_throw_ia(x, y, index_min_max(2, 3), index_min(10));
test_throw_ia(x, y, index_min_max(1, 3), index_min(2));
}

Expand Down
8 changes: 4 additions & 4 deletions src/test/unit/model/indexing/assign_varmat_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,7 @@ void test_min_vec() {
check_adjs(check_i, x, "lhs");
check_adjs([](int /* i */) { return true; }, y, "rhs");
test_throw_out_of_range(x, y, index_min(0));
test_throw_out_of_range(x, y, index_min(6));
test_throw_invalid_arg(x, y, index_min(6));
test_throw_invalid_arg(x, conditionally_generate_linear_var_vector<Vec>(4),
index_min(3));
test_throw_invalid_arg(x, conditionally_generate_linear_var_vector<Vec>(2),
Expand Down Expand Up @@ -1104,7 +1104,7 @@ void min_matrix_test() {
check_adjs(check_i_x, check_all, x, "lhs", 0);
check_adjs(check_all, y, "rhs", 1.0);
test_throw_out_of_range(x, y, index_min(0));
test_throw_out_of_range(x, y, index_min(4));
test_throw_invalid_arg(x, y, index_min(4));
test_throw_invalid_arg(x, y, index_min(1));
var_value<MatrixXd> z(MatrixXd::Ones(1, 2));
test_throw_invalid_arg(x, z, index_min(2));
Expand Down Expand Up @@ -1137,7 +1137,7 @@ void minmax_min_matrix_test() {
test_throw_out_of_range(x, y, index_min_max(0, 3), index_min(2));
test_throw_out_of_range(x, y, index_min_max(2, 4), index_min(2));
test_throw_out_of_range(x, y, index_min_max(2, 3), index_min(0));
test_throw_out_of_range(x, y, index_min_max(2, 3), index_min(5));
test_throw_invalid_arg(x, y, index_min_max(2, 3), index_min(5));
test_throw_invalid_arg(x, conditionally_generate_linear_var_matrix(1, 3, 10),
index_min_max(2, 3), index_min(2));
test_throw_invalid_arg(x, conditionally_generate_linear_var_matrix(2, 5, 10),
Expand Down Expand Up @@ -1202,7 +1202,7 @@ void min_max_matrix_test() {
auto check_all = [](int /* i*/) { return true; };
check_adjs(check_all, check_all, y, "rhs");
test_throw_out_of_range(x, y, index_min(0), index_max(2));
test_throw_out_of_range(x, y, index_min(5), index_max(2));
test_throw_invalid_arg(x, y, index_min(5), index_max(2));
test_throw_invalid_arg(x, y, index_min(2), index_max(0));
test_throw_out_of_range(x, y, index_min(2), index_max(5));
test_throw_invalid_arg(x, y, index_min(2), index_max(1));
Expand Down
180 changes: 180 additions & 0 deletions src/test/unit/model/indexing/empty_range_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,180 @@
#include <stan/model/indexing.hpp>
#include <stan/math/rev/fun/sum.hpp>
#include <stan/math/rev/fun/value_of.hpp>
#include <gtest/gtest.h>
#include <limits>
#include <vector>

using stan::model::assign;
using stan::model::index_max;
using stan::model::index_min;
using stan::model::index_multi;
using stan::model::index_omni;
using stan::model::index_uni;
using stan::model::rvalue;

namespace {

// An empty read contributes no adjoints; an empty assignment leaves the
// original values and their derivatives unchanged.
template <typename T, typename... Idxs>
void check_empty_slice(T& x, int rows, int cols, const Idxs&... idxs) {
auto selected = rvalue(x, "x", idxs...);
EXPECT_EQ(rows, selected.rows());
EXPECT_EQ(cols, selected.cols());
if constexpr (stan::is_var<stan::scalar_type_t<T>>::value) {
stan::math::set_zero_all_adjoints();
stan::math::sum(selected).grad();
for (int i = 0; i < x.size(); ++i) {
if constexpr (stan::is_var_matrix<T>::value) {
EXPECT_EQ(0, x.adj().coeff(i));
} else {
EXPECT_EQ(0, x.coeff(i).adj());
}
}
}
using plain_t = stan::plain_type_t<decltype(selected)>;
plain_t empty(selected);
EXPECT_NO_THROW(assign(x, empty, "x", idxs...));
EXPECT_TRUE(stan::math::value_of(x).isOnes());
if constexpr (stan::is_var<stan::scalar_type_t<T>>::value) {
stan::math::set_zero_all_adjoints();
stan::math::sum(x).grad();
for (int i = 0; i < x.size(); ++i) {
if constexpr (stan::is_var_matrix<T>::value) {
EXPECT_EQ(1, x.adj().coeff(i));
} else {
EXPECT_EQ(1, x.coeff(i).adj());
}
}
}
}

template <typename T>
void check_vector_empty_ranges() {
using values_t
= Eigen::Matrix<double, T::RowsAtCompileTime, T::ColsAtCompileTime>;
for (int size : {0, 3}) {
T x(values_t::Ones(size));
const int rows = T::RowsAtCompileTime == 1 ? 1 : 0;
const int cols = T::RowsAtCompileTime == 1 ? 0 : 1;
for (int min : {size + 1, size + 4, std::numeric_limits<int>::max()}) {
check_empty_slice(x, rows, cols, index_min(min));
EXPECT_THROW(assign(x, values_t::Ones(1), "x", index_min(min)),
std::invalid_argument);
}
for (int max : {0, -3, std::numeric_limits<int>::min()}) {
check_empty_slice(x, rows, cols, index_max(max));
EXPECT_THROW(assign(x, values_t::Ones(1), "x", index_max(max)),
std::invalid_argument);
}
EXPECT_THROW(rvalue(x, "x", index_min(0)), std::out_of_range);
EXPECT_THROW(rvalue(x, "x", index_max(size + 1)), std::out_of_range);
}
}

template <typename T>
void check_matrix_empty_ranges() {
for (int rows : {0, 3}) {
for (int cols : {0, 4}) {
T x(Eigen::MatrixXd::Ones(rows, cols));
for (int min : {rows + 1, rows + 4, std::numeric_limits<int>::max()}) {
check_empty_slice(x, 0, cols, index_min(min));
check_empty_slice(x, 0, cols, index_min(min), index_omni());
EXPECT_THROW(assign(x, Eigen::MatrixXd(1, cols), "x", index_min(min)),
std::invalid_argument);
EXPECT_THROW(
assign(x, Eigen::MatrixXd(0, cols + 1), "x", index_min(min)),
std::invalid_argument);
if (cols > 0) {
check_empty_slice(x, 0, 1, index_min(min), index_uni(1));
check_empty_slice(x, 0, 2, index_min(min),
index_multi(std::vector<int>{1, 2}));
}
}
for (int min : {cols + 1, cols + 4, std::numeric_limits<int>::max()}) {
check_empty_slice(x, rows, 0, index_omni(), index_min(min));
check_empty_slice(x, 0, 0, index_min(rows + 4), index_min(min));
check_empty_slice(x, 0, 0, index_max(0), index_min(min));
EXPECT_THROW(assign(x, Eigen::MatrixXd(rows, 1), "x", index_omni(),
index_min(min)),
std::invalid_argument);
if (rows > 0) {
check_empty_slice(x, 1, 0, index_uni(1), index_min(min));
check_empty_slice(x, 2, 0, index_multi(std::vector<int>{1, 2}),
index_min(min));
}
}
for (int max : {0, -3, std::numeric_limits<int>::min()}) {
check_empty_slice(x, 0, cols, index_max(max));
check_empty_slice(x, rows, 0, index_omni(), index_max(max));
check_empty_slice(x, 0, 0, index_min(rows + 4), index_max(max));
EXPECT_THROW(assign(x, Eigen::MatrixXd(1, cols), "x", index_max(max)),
std::invalid_argument);
EXPECT_THROW(assign(x, Eigen::MatrixXd(rows, 1), "x", index_omni(),
index_max(max)),
std::invalid_argument);
}
EXPECT_THROW(rvalue(x, "x", index_min(0)), std::out_of_range);
EXPECT_THROW(rvalue(x, "x", index_omni(), index_min(0)),
std::out_of_range);
EXPECT_THROW(rvalue(x, "x", index_max(rows + 1)), std::out_of_range);
EXPECT_THROW(rvalue(x, "x", index_omni(), index_max(cols + 1)),
std::out_of_range);
}
}
}

} // namespace

TEST(ModelIndexingEmptyRange, eigen) {
check_vector_empty_ranges<Eigen::VectorXd>();
check_vector_empty_ranges<Eigen::RowVectorXd>();
check_matrix_empty_ranges<Eigen::MatrixXd>();
}

TEST(ModelIndexingEmptyRange, eigenVar) {
stan::math::nested_rev_autodiff nested;
check_vector_empty_ranges<Eigen::Matrix<stan::math::var, -1, 1>>();
check_vector_empty_ranges<Eigen::Matrix<stan::math::var, 1, -1>>();
check_matrix_empty_ranges<Eigen::Matrix<stan::math::var, -1, -1>>();
}

TEST(ModelIndexingEmptyRange, varmat) {
stan::math::nested_rev_autodiff nested;
check_vector_empty_ranges<stan::math::var_value<Eigen::VectorXd>>();
check_vector_empty_ranges<stan::math::var_value<Eigen::RowVectorXd>>();
check_matrix_empty_ranges<stan::math::var_value<Eigen::MatrixXd>>();
}

TEST(ModelIndexingEmptyRange, arrays) {
for (int size : {0, 3}) {
std::vector<double> x(size, 1);
std::vector<std::vector<double>> xx(size, x);
const auto check = [&](auto idx) {
EXPECT_TRUE(rvalue(x, "x", idx).empty());
EXPECT_TRUE(rvalue(xx, "xx", idx, index_omni()).empty());
EXPECT_NO_THROW(assign(x, std::vector<double>{}, "x", idx));
EXPECT_NO_THROW(assign(xx, std::vector<std::vector<double>>{}, "xx", idx,
index_omni()));
EXPECT_THROW(assign(x, std::vector<double>{2}, "x", idx),
std::invalid_argument);
if (size > 0) {
auto inner_empty = rvalue(xx, "xx", index_omni(), idx);
EXPECT_EQ(size, inner_empty.size());
for (const auto& inner : inner_empty) {
EXPECT_TRUE(inner.empty());
}
EXPECT_NO_THROW(assign(xx, inner_empty, "xx", index_omni(), idx));
}
EXPECT_EQ(std::vector<double>(size, 1), x);
EXPECT_EQ(std::vector<std::vector<double>>(size, x), xx);
};
for (int min : {size + 1, size + 4, std::numeric_limits<int>::max()}) {
check(index_min(min));
}
for (int max : {0, -3, std::numeric_limits<int>::min()}) {
check(index_max(max));
}
}
}
Loading
Loading