Skip to content

add overloads for scalar complex and eigen matrix multiplication #2680

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
SteveBronder merged 4 commits into from
Mar 2, 2022
Merged

add overloads for scalar complex and eigen matrix multiplication #2680

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
de84a2a
add overloads for scalar complex and eigen matrix multiplication
SteveBronder Feb 28, 2022
4217894
update headers
SteveBronder Feb 28, 2022
687fea2
update docs
SteveBronder Feb 28, 2022
ad4ecc5
add complex ad tests
SteveBronder Mar 1, 2022
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
8 changes: 5 additions & 3 deletions stan/math/prim/core/operator_multiplication.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
#define STAN_MATH_PRIM_CORE_OPERATOR_MULTIPLICATION_HPP

#include <stan/math/prim/meta/return_type.hpp>
#include <stan/math/prim/meta/is_complex.hpp>
#include <stan/math/prim/meta/is_stan_scalar.hpp>
#include <complex>

namespace stan {
Expand All @@ -18,7 +20,7 @@ namespace internal {
* @param[in] rhs second argument
* @return product of the arguments
*/
template <typename U, typename V>
template <typename U, typename V, require_any_complex_t<U, V>* = nullptr>
inline complex_return_t<U, V> complex_multiply(const U& lhs, const V& rhs) {
complex_return_t<U, V> y(lhs);
y *= rhs;
Expand Down Expand Up @@ -50,7 +52,7 @@ inline complex_return_t<U, V> operator*(const std::complex<U>& x,
* @param y second argument
* @return product of the arguments
*/
template <typename U, typename V>
template <typename U, typename V, require_stan_scalar_t<V>* = nullptr>
inline complex_return_t<U, V> operator*(const std::complex<U>& x, const V& y) {
return internal::complex_multiply(x, y);
}
Expand All @@ -64,7 +66,7 @@ inline complex_return_t<U, V> operator*(const std::complex<U>& x, const V& y) {
* @param y second argument
* @return product of the arguments
*/
template <typename U, typename V>
template <typename U, typename V, require_stan_scalar_t<U>* = nullptr>
inline complex_return_t<U, V> operator*(const U& x, const std::complex<V>& y) {
return internal::complex_multiply(x, y);
}
Expand Down
42 changes: 40 additions & 2 deletions stan/math/prim/fun/multiply.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,48 @@ namespace math {
*/
template <typename Mat, typename Scal, require_stan_scalar_t<Scal>* = nullptr,
require_eigen_t<Mat>* = nullptr,
require_all_not_st_var<Scal, Mat>* = nullptr>
require_all_not_st_var<Scal, Mat>* = nullptr,
require_all_not_complex_t<Scal, value_type_t<Mat>>* = nullptr>
inline auto multiply(const Mat& m, Scal c) {
return c * m;
}

/**
* Return the product of the specified matrix and scalar, one of which must have a complex value type.
* The return type will be an expression template denoting a complex matrix.
*
* @tparam Mat type of matrix argument
* @tparam Scal type of scalar argument
* @param m matrix argument
* @param c scalar argument
* @return expression template evaluating to the product of the matrix and scalar arguments
*/
template <typename Mat, typename Scal,
require_any_complex_t<value_type_t<Mat>, Scal>* = nullptr,
require_eigen_t<Mat> * = nullptr,
require_not_eigen_t<Scal>* = nullptr>
inline auto multiply(const Mat& m, Scal c) {
return m * c;
}

/**
* Return the product of the specified matrix and scalar, one of which must have a complex value type.
* The return type will be an expression template denoting a complex matrix.
*
* @tparam Scal type of scalar argument
* @tparam Mat type of matrix argument
* @param c scalar argument
* @param m matrix argument
* @return expression template evaluating to the product of the matrix and scalar arguments
*/
template <typename Mat, typename Scal,
require_any_complex_t<value_type_t<Mat>, Scal>* = nullptr,
require_eigen_t<Mat> * = nullptr,
require_not_eigen_t<Scal>* = nullptr>
inline auto multiply(const Scal& m, const Mat& c) {
return m * c;
}

/**
* Return specified scalar multiplied by specified matrix.
*
Expand All @@ -39,7 +76,8 @@ inline auto multiply(const Mat& m, Scal c) {
*/
template <typename Scal, typename Mat, require_stan_scalar_t<Scal>* = nullptr,
require_eigen_t<Mat>* = nullptr,
require_all_not_st_var<Scal, Mat>* = nullptr>
require_all_not_st_var<Scal, Mat>* = nullptr,
require_all_not_complex_t<Scal, value_type_t<Mat>>* = nullptr>
inline auto multiply(Scal c, const Mat& m) {
return c * m;
}
Expand Down
2 changes: 2 additions & 0 deletions stan/math/rev/fun/multiply.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,6 +189,8 @@ inline auto multiply(const T1& A, const T2& B) {
template <typename T1, typename T2, require_matrix_t<T1>* = nullptr,
require_not_matrix_t<T2>* = nullptr,
require_any_st_var<T1, T2>* = nullptr,
require_not_complex_t<value_type_t<T1>>* = nullptr,
require_not_complex_t<value_type_t<T2>>* = nullptr,
require_not_row_and_col_vector_t<T1, T2>* = nullptr>
inline auto multiply(const T1& A, const T2& B) {
return multiply(B, A);
Expand Down
50 changes: 50 additions & 0 deletions test/unit/math/mix/fun/multiply_complex_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
#include <test/unit/math/test_ad.hpp>
#include <test/unit/math/mix/fun/multiply_util.hpp>

TEST(mathMix, complexMultiply) {
auto f
= [](const auto& x, const auto& y) { return stan::math::multiply(x, y); };
double d_scalar = 1.0;
std::complex<double> cd_scalar = 1.0;
Eigen::Matrix<double, -1, -1> d_mat(2, 2);
d_mat << 1, 2, 3, 4;
Eigen::Matrix<std::complex<double>, -1, -1> cd_mat(2, 2);
cd_mat << 1, 2, 3, 4;

stan::test::expect_ad(f, d_scalar, cd_mat);
stan::test::expect_ad(f, cd_scalar, d_mat);
stan::test::expect_ad(f, cd_scalar, cd_mat);
stan::test::expect_ad(f, cd_mat, d_scalar);
stan::test::expect_ad(f, d_mat, cd_scalar);
stan::test::expect_ad(f, cd_mat, cd_scalar);
stan::test::expect_ad(f, d_mat, cd_mat);
stan::test::expect_ad(f, cd_mat, d_mat);
stan::test::expect_ad(f, cd_mat, cd_mat);

Eigen::Matrix<double, -1, 1> d_vec(2);
d_vec << 1, 2;
Eigen::Matrix<std::complex<double>, -1, 1> cd_vec(2);
cd_vec << 1, 2;

stan::test::expect_ad(f, d_mat, cd_vec);
stan::test::expect_ad(f, cd_mat, d_vec);
stan::test::expect_ad(f, cd_mat, cd_vec);

Eigen::Matrix<double, 1, -1> d_rowvec(2);
d_rowvec << 1, 2;
Eigen::Matrix<std::complex<double>, 1, -1> cd_rowvec(2);
cd_rowvec << 1, 2;

stan::test::expect_ad(f, d_vec, cd_rowvec);
stan::test::expect_ad(f, cd_vec, d_rowvec);
stan::test::expect_ad(f, cd_vec, cd_rowvec);
stan::test::expect_ad(f, d_rowvec, cd_vec);
stan::test::expect_ad(f, cd_rowvec, d_vec);
stan::test::expect_ad(f, cd_rowvec, cd_vec);
stan::test::expect_ad(f, d_rowvec, cd_mat);
stan::test::expect_ad(f, cd_rowvec, d_mat);
stan::test::expect_ad(f, cd_rowvec, cd_mat);
stan::test::expect_ad(f, d_mat, cd_vec);
stan::test::expect_ad(f, cd_mat, d_vec);
stan::test::expect_ad(f, cd_mat, cd_vec);
}
44 changes: 44 additions & 0 deletions test/unit/math/mix/fun/multiply_util.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,10 @@
template <typename T>
void instantiate_multiply() {
using stan::math::multiply;
T v_scalar = 1.0;
double d_scalar = 1.0;
std::complex<double> cd_scalar = 1.0;
std::complex<T> cv_scalar = 1.0;
Eigen::Matrix<double, -1, -1> d_mat(2, 2);
d_mat << 1, 2, 3, 4;
Eigen::Matrix<T, -1, -1> v_mat(2, 2);
Expand All @@ -14,6 +18,46 @@ void instantiate_multiply() {
Eigen::Matrix<std::complex<T>, -1, -1> cv_mat(2, 2);
cv_mat << 1, 2, 3, 4;

auto d_scalar_d_mat = stan::math::eval(multiply(d_scalar, d_mat));

@bob-carpenter bob-carpenter Feb 28, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[question]
Does this run autodiff and value tests? If not, we should also have those tests.

@SteveBronder SteveBronder Feb 28, 2022

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It doesn't because we currently do not support complex types in the ad testing framework but in a seperate PR I can add that. I think we just need to add a compile time flag to the ad testing scheme where if it's true then we also run the function with complex types

auto d_scalar_v_mat = stan::math::eval(multiply(d_scalar, v_mat));
auto d_scalar_cd_mat = stan::math::eval(multiply(d_scalar, cd_mat));
auto d_scalar_cv_mat = stan::math::eval(multiply(d_scalar, cv_mat));

auto v_scalar_d_mat = stan::math::eval(multiply(v_scalar, d_mat));
auto v_scalar_v_mat = stan::math::eval(multiply(v_scalar, v_mat));
auto v_scalar_cd_mat = stan::math::eval(multiply(v_scalar, cd_mat));
auto v_scalar_cv_mat = stan::math::eval(multiply(v_scalar, cv_mat));

auto cd_scalar_d_mat = stan::math::eval(multiply(cd_scalar, d_mat));
auto cd_scalar_v_mat = stan::math::eval(multiply(cd_scalar, v_mat));
auto cd_scalar_cd_mat = stan::math::eval(multiply(cd_scalar, cd_mat));
auto cd_scalar_cv_mat = stan::math::eval(multiply(cd_scalar, cv_mat));

auto cv_scalar_d_mat = stan::math::eval(multiply(cv_scalar, d_mat));
auto cv_scalar_v_mat = stan::math::eval(multiply(cv_scalar, v_mat));
auto cv_scalar_cd_mat = stan::math::eval(multiply(cv_scalar, cd_mat));
auto cv_scalar_cv_mat = stan::math::eval(multiply(cv_scalar, cv_mat));

auto d_mat_d_scalar = stan::math::eval(multiply(d_mat, d_scalar));
auto v_mat_d_scalar = stan::math::eval(multiply(v_mat, d_scalar));
auto cd_mat_d_scalar = stan::math::eval(multiply(cd_mat, d_scalar));
auto cv_mat_d_scalar = stan::math::eval(multiply(cv_mat, d_scalar));

auto d_mat_v_scalar = stan::math::eval(multiply(d_mat, v_scalar));
auto v_mat_v_scalar = stan::math::eval(multiply(v_mat, v_scalar));
auto cd_mat_v_scalar = stan::math::eval(multiply(cd_mat, v_scalar));
auto cv_mat_v_scalar = stan::math::eval(multiply(cv_mat, v_scalar));

auto d_mat_cd_scalar = stan::math::eval(multiply(d_mat, cd_scalar));
auto v_mat_cd_scalar = stan::math::eval(multiply(v_mat, cd_scalar));
auto cd_mat_cd_scalar = stan::math::eval(multiply(cd_mat, cd_scalar));
auto cv_mat_cd_scalar = stan::math::eval(multiply(cv_mat, cd_scalar));

auto d_mat_cv_scalar = stan::math::eval(multiply(d_mat, cv_scalar));
auto v_mat_cv_scalar = stan::math::eval(multiply(v_mat, cv_scalar));
auto cd_mat_cv_scalar = stan::math::eval(multiply(cd_mat, cv_scalar));
auto cv_mat_cv_scalar = stan::math::eval(multiply(cv_mat, cv_scalar));

auto d_d_mat = stan::math::eval(multiply(d_mat, d_mat));
auto d_v_mat = stan::math::eval(multiply(d_mat, v_mat));
auto d_cd_mat = stan::math::eval(multiply(d_mat, cd_mat));
Expand Down