See More

// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*- // [[Rcpp::depends(BH)]] // [[Rcpp::depends(RcppEigen)]] // [[Rcpp::depends(StanHeaders)]] #include #include #include // [[Rcpp::plugins(cpp14)]] using Eigen::Map; using Eigen::Dynamic; using Eigen::Matrix; //The following template needs to be defined to do AD struct scoretemp { //Following inputs are all constant w.r.t. differentiation const Matrix S_; // Linear constraints matrix const Matrix y_; // Realisation const Matrix x_; // Draws from probabilistic forecast const Matrix xs_; // Draws from probabilistic forecast int score_;//Score to be used double alpha_; //Additional parameters scoretemp(const Matrix& S, const Matrix& y, const Matrix& x, const Matrix& xs, int score, double alpha) : S_(S),y_(y),x_(x),xs_(xs),score_(score),alpha_(alpha) { } template T operator()(const Matrix& Gvec) const{ int n = S_.rows(); // Number of series int m = S_.cols(); // Number of base series int Q = x_.cols(); // Number of draws from probabilistic forecast Matrix g1 = Gvec.topRows(m); //First m parameters for translation Matrix g2 = Gvec.bottomRows(n*m); //Remaining parameters for matrix Matrix G = to_matrix(g2 , m , n); //Put into matrix Matrix SG = stan::math::multiply(S_,G); // S times G (needs stan::math::multiply) Matrix xr = stan::math::multiply(SG,x_); // Draw from reconciled forecast (needs stan::math::multiply) Matrix xrs = stan::math::multiply(SG,xs_); // Draw from reconciled forecast (needs stan::math::multiply) Matrix d = stan::math::multiply(S_,g1); // Translation T out; // output if (score_==1){ T term1 = 0; // Intialise terms for cumulative sum T term2 = 0; // Intialise terms for sumulative sum Matrix yd = y_-d; // y-d Matrix dif1; // Initialise vectors Matrix dif2; // Initialise vectors for (int q=0; q> S(Rcpp::as> >(Sin)); Map> y(Rcpp::as> >(yin)); Map> x(Rcpp::as> >(xin)); Map> xs(Rcpp::as> >(xsin)); Map> G(Rcpp::as> >(Gin)); scoretemp f(S,y,x,xs,scorein,alphain); // Instantiate template double fx; // Value of function Matrix grad_fx; // Initialise gradient stan::math::gradient(f, G, fx, grad_fx); //Automatic differentiation // Convert back to list return Rcpp::List::create(Rcpp::Named("grad") = Rcpp::wrap(grad_fx), Rcpp::Named("val") = fx); }