Hi,
I have been working on learning how to implement my own distribution that uses some optimised c++ code I already have. I have worked through the process of implementing a copy of the normal distribution successfully and have looked at the multi normal cholesky for some guidance on other things. I have gotten to the point where I can call my optimised c++ function (which takes pointers to doubles and returns the log likelihood and gradients via pointers as well). However, when I try to assign the computed gradients to ops_partials.edge1_.partials_ I get an error about undefined reference to empty broadcast array. If I don’t set anything to ops_partials then it compiles fine, so i think its something wrong with how I am setting those values.
Here is the code for my new custom distribution, I have actually commented out the bits that would
#ifndef STAN_MATH_PRIM_PROB_MATRIX_NORMAL_CHOLESKY_LPDF_HPP
#define STAN_MATH_PRIM_PROB_MATRIX_NORMAL_CHOLESKY_LPDF_HPP
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/as_value_column_vector_or_scalar.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/dot_self.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/max_size_mvt.hpp>
#include <stan/math/prim/fun/mdivide_left_tri.hpp>
#include <stan/math/prim/fun/mdivide_right_tri.hpp>
#include <stan/math/prim/fun/size_mvt.hpp>
#include <stan/math/prim/fun/sum.hpp>
#include <stan/math/prim/fun/transpose.hpp>
#include <stan/math/prim/fun/vector_seq_view.hpp>
#include <stan/math/prim/functor/operands_and_partials.hpp>
#include<lib/eigen_3.3.9/Eigen/src/misc/blas.h>
#include "ll_dll_normal.hpp"
namespace stan {
namespace math {
/** \ingroup multivar_dists
* The log of the matrix normal density for the given y, mu, and LSigma
* where L is the lower cholesky factor of a covariance matrix such that LSigma*LSigma' = Sigma
*
* @tparam T_y type of scalar
* @tparam T_mu type of location
* @tparam T_LSigma type of LSigma
*
* @param y An mxn matrix.
* @param mu The mean matrix.
* @param LSigma The nxn the lower cholesky factor of the covariance matrix Sigma
* @return The log of the matrix normal density.
* @throw std::domain_error if Sigma is not square
*/
template <bool propto, typename T_y, typename T_mu, typename T_LSigma,
require_all_matrix_t<T_y, T_mu, T_LSigma >* = nullptr>
return_type_t<T_y, T_mu, T_LSigma> matrix_normal_cholesky_lpdf(
const T_y& y, const T_mu& mu, const T_LSigma& LSigma) {
static const char* function = "matrix_normal_cholesky_lpdf";
check_finite(function, "LSigma", LSigma);
check_size_match(function, "Rows of random variable", y.rows(),
"Rows of location parameter", mu.rows());
check_size_match(function, "Columns of random variable", y.cols(),
"Columns of location parameter", mu.cols());
check_size_match(function, "Rows of random variable", y.rows(),
"Rows of LSigma", LSigma.rows());
check_finite(function, "Location parameter", mu);
check_finite(function, "Random variable", y);
using T_covar_elem = typename scalar_type<T_LSigma>::type;
using T_return = return_type_t<T_y, T_mu, T_LSigma>;
using T_partials_return = partials_return_t<T_y, T_mu, T_LSigma>;
using matrix_partials_t
= Eigen::Matrix<T_partials_return, Eigen::Dynamic, Eigen::Dynamic>;
using T_y_ref = ref_type_t<T_y>;
using T_mu_ref = ref_type_t<T_mu>;
using T_L_ref = ref_type_t<T_LSigma>;
T_y_ref y_ref = y;
T_mu_ref mu_ref = mu;
T_L_ref L_ref = LSigma;
decltype(auto) y_val = to_ref(value_of(y_ref));
decltype(auto) mu_val = to_ref(value_of(mu_ref));
decltype(auto) L_val = to_ref(value_of(L_ref));
// prepare data for external call
int p = y.rows();
int n = y.cols();
double *mu_data = const_cast<double*>(mu_val.data());
double *y_data = const_cast<double*>(y_val.data());
double *L_data = const_cast<double*>(L_val.data());
// get only the triangular bits of L_data that we want
int ntheta = p*(p+1)/2;
double * theta = (double *)malloc(ntheta*sizeof(double));
int count = 0;
for (size_t i = 0; i < p; i++){
for (size_t j = 0; j < i+1; j++){
theta[count] = L_data[i,j];
count++;
}
}
double * dLL = (double *)malloc(p*n+p*n+ntheta*sizeof(double));
double LL = 0;
blas_maths::ll_and_dll(p, n, y_data, mu_data, theta, &LL, dLL);
int ione = 1;
int itmp = n*p;
matrix_partials_t dLLdy(p, n);
// dLLdy = matrix_partials_t::Zero(p, n);
dcopy_(&itmp,dLL,&ione,dLLdy.data(),&ione);
matrix_partials_t dLLdmu(p, n);
dLLdmu = matrix_partials_t::Zero(p, n);
// dcopy_(&itmp,dLL+n*p,&ione,dLLdmu.data(),&ione);
matrix_partials_t dLLdLSigma(p, p);
dLLdLSigma = matrix_partials_t::Zero(p, p); // zero this as we will only then set lower triangle
dcopy_(&itmp,dLL+n*p,&ione,dLLdmu.data(),&ione);
count = 0;
for (size_t i = 0; i < p; i++){
for (size_t j = 0; j < i+1; j++){
dLLdLSigma[i, j] = dLL[n*p+n*p+count];
count++;
}
}
operands_and_partials<T_y_ref, T_mu_ref, T_L_ref> ops_partials(y_ref, mu_ref, L_ref);
T_partials_return logp(0);
matrix_partials_t inv_L_val = mdivide_left_tri<Eigen::Lower>(value_of(L_ref));
logp = LL;
ops_partials.edge1_.partials_ = dLLdy;
ops_partials.edge2_.partials_ = dLLdmu;
ops_partials.edge3_.partials_ = dLLdLSigma;
return ops_partials.build(logp);
}
template <typename T_y, typename T_mu, typename T_LSigma,
require_all_matrix_t<T_y, T_mu, T_LSigma>* = nullptr>
return_type_t<T_y, T_mu, T_LSigma> matrix_normal_cholesky_lpdf(
const T_y& y, const T_mu& mu, const T_LSigma& LSigma) {
return matrix_normal_cholesky_lpdf<false>(y, mu, LSigma);
}
} // namespace math
} // namespace stan
#endif
The error I am getting is
/usr/bin/ld: ../models/lpv_v01_MNC.o: in function `stan::return_type<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >::type stan::math::matrix_normal_cholesky_lpdf<false, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, (void*)0>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&)':
lpv_v01_MNC.hpp:(.text._ZN4stan4math27matrix_normal_cholesky_lpdfILb0EN5Eigen6MatrixIdLin1ELin1ELi0ELin1ELin1EEES4_S4_LPv0EEENS_11return_typeIJT0_T1_T2_EE4typeERKS7_RKS8_RKS9_[_ZN4stan4math27matrix_normal_cholesky_lpdfILb0EN5Eigen6MatrixIdLin1ELin1ELi0ELin1ELin1EEES4_S4_LPv0EEENS_11return_typeIJT0_T1_T2_EE4typeERKS7_RKS8_RKS9_]+0x868): undefined reference to `stan::math::internal::empty_broadcast_array<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, void>::operator+=(Eigen::Matrix<double, -1, -1, 0, -1, -1>)'
/usr/bin/ld: ../models/lpv_v01_MNC.o: in function `stan::return_type<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >::type stan::math::matrix_normal_cholesky_lpdf<true, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, (void*)0>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&)':
lpv_v01_MNC.hpp:(.text._ZN4stan4math27matrix_normal_cholesky_lpdfILb1EN5Eigen6MatrixIdLin1ELin1ELi0ELin1ELin1EEES4_S4_LPv0EEENS_11return_typeIJT0_T1_T2_EE4typeERKS7_RKS8_RKS9_[_ZN4stan4math27matrix_normal_cholesky_lpdfILb1EN5Eigen6MatrixIdLin1ELin1ELi0ELin1ELin1EEES4_S4_LPv0EEENS_11return_typeIJT0_T1_T2_EE4typeERKS7_RKS8_RKS9_]+0x868): undefined reference to `stan::math::internal::empty_broadcast_array<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, void>::operator+=(Eigen::Matrix<double, -1, -1, 0, -1, -1>)'
collect2: error: ld returned 1 exit status
Here is the external c++ code in case that helps
#ifndef LL_DLL_NORMAL_HPP
#define LL_DLL_NORMAL_HPP
#include <string.h>
#include<lib/eigen_3.3.9/Eigen/src/misc/blas.h>
#include<lib/eigen_3.3.9/Eigen/src/misc/lapack.h>
#include "math.h"
namespace blas_maths {
void ll_and_dll(int p,
int N,
double *y,
double *yh,
double *theta,
double *LL,
double *dLL)
{
double *pe, *Pi, *R, *S, dtmp;
int k, i, j, itmp;
// Extract Gamma and Pi from theta
Pi = theta;
int ione = 1;
int izero = 0;
double done = 1.0;
double dzero = 0.0;
double dmone = -1.0;
char LOW = 'L';
char UPP = 'U';
char TRN = 'T';
char NO = 'N';
char LEFT = 'L';
char RIGHT = 'R';
//Create memory and zero it
pe = (double *)malloc(p*N*sizeof(double));
R = (double *)malloc(p*p*sizeof(double));
S = (double *)malloc(p*p*sizeof(double));
dtmp = 0.0;
itmp = p*p;
dcopy_(&itmp,&dtmp,&izero,R,&ione);
/* Clear R and set it up */
k = 0;
for(i=0;i<p;i++){
for(j=0;j<=i;j++){
R[j + i*p] = Pi[k];
k++;
}
}
//Form error: pe=y-yh
itmp = p*N;
dtmp = 0.0;
dcopy_(&itmp,y,&ione,pe,&ione);
daxpy_(&itmp,&dmone,yh,&ione,pe,&ione);
//Scale error by inverse transpose Cholesky factor: pe = R' \ pe
dtrsm_(&LEFT,&UPP,&TRN,&NO,&p,&N,&done,R,&p,pe,&p);
//Form S = (R\pe)*(R\pe)'
dsyrk_(&UPP, &NO, &p, &N, &done, pe, &p, &dzero, S, &p);
for(i=0;i<p;i++){
for(j=0;j<i;j++){
S[i+j*p] = S[j+i*p];
}
}
//Compute LL
LL[0] = 0.0;
for(i=0;i<p;i++){
LL[0] = LL[0] - (0.5*N)*log(R[i+i*p]*R[i+i*p]) - 0.5*S[i+i*p];
}
//Compute gradient
double *dLLdy, *dLLdyh, *dLLdR;
dLLdy = dLL;
i = p*N;
dLLdyh = &dLL[i];
i = i+p*N;
dLLdR = &dLL[i];
//pe = R \ (R' \ (y-yh))
dtrsm_(&LEFT,&UPP,&NO,&NO,&p,&N,&done,R,&p,pe,&p);
itmp = p*N;
dcopy_(&itmp,&dzero,&izero,dLLdy,&ione);
daxpy_(&itmp,&dmone,pe,&ione,dLLdy,&ione);
dcopy_(&itmp,pe,&ione,dLLdyh,&ione);
// 2*vech(diag(Z.Ny./diag(R)) - (S/R')
dtrsm_(&RIGHT,&UPP,&TRN,&NO,&p,&p,&done,R,&p,S,&p);
k=0;
for(i=0;i<p;i++){
S[i+i*p] = -N/R[i+i*p] + S[i+i*p];
for(j=0;j<=i;j++){
dLLdR[k] = S[j+i*p];
k++;
}
}
exit_stub:
free(pe);
free(R);
free(S);
}
}
#endif
//void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
// double *y, *yh, *theta, *LL, *dLL, *yout;
// int p, N, ntheta;
//
//
// // n, m, p, na, nb, N, smoothing, y, u, a, b, theta, P1, X1
// y = mxGetPr(prhs[0]);
// p = mxGetM(prhs[0]);
// N = mxGetN(prhs[0]);
// yh = mxGetPr(prhs[1]);
// theta = mxGetPr(prhs[2]);
// ntheta = mxGetM(prhs[2]);
//
// plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
// LL = mxGetPr(plhs[0]);
// plhs[1] = mxCreateDoubleMatrix(p*N+p*N+ntheta,1,mxREAL);
// dLL = mxGetPr(plhs[1]);
//
// ll_and_dll(p,N,y,yh,theta,LL,dLL);
//}
Thanks in advance for any help. I am feeling very stuck on this. Also, if there is any advice on better ways to access the data as point to double and to set it back into eigen matrices then that would be great.