##### https://github.com/cran/rstpm2
Tip revision: 11f1ef1
multistate.cpp
``````#include <RcppArmadillo.h>
#include <errno.h>
#include <math.h>
// #include <boost/numeric/odeint.hpp>
#include <vector>
#include <algorithm>

/*
*	Natural Splines
*	---------------
*	Here the end-conditions are determined by setting the second
*	derivative of the spline at the end-points to equal to zero.
*
*	There are n-2 unknowns (y[i]'' at x[2], ..., x[n-1]) and n-2
*	equations to determine them.  Either Choleski or Gaussian
*	elimination could be used.
*/
static void
natural_spline(int n, double *x, double *y, double *b, double *c, double *d)
{
int nm1, i;
double t;
x--; y--; b--; c--; d--;
if(n < 2) {
errno = EDOM;
return;
}
if(n < 3) {
t = (y[2] - y[1]);
b[1] = t / (x[2]-x[1]);
b[2] = b[1];
c[1] = c[2] = d[1] = d[2] = 0.0;
return;
}
nm1 = n - 1;
/* Set up the tridiagonal system */
/* b = diagonal, d = offdiagonal, c = right hand side */
d[1] = x[2] - x[1];
c[2] = (y[2] - y[1])/d[1];
for( i=2 ; i<n ; i++) {
d[i] = x[i+1] - x[i];
b[i] = 2.0 * (d[i-1] + d[i]);
c[i+1] = (y[i+1] - y[i])/d[i];
c[i] = c[i+1] - c[i];
}
/* Gaussian elimination */
for(i=3 ; i<n ; i++) {
t = d[i-1]/b[i-1];
b[i] = b[i] - t*d[i-1];
c[i] = c[i] - t*c[i-1];
}
/* Backward substitution */
c[nm1] = c[nm1]/b[nm1];
for(i=n-2 ; i>1 ; i--)
c[i] = (c[i]-d[i]*c[i+1])/b[i];
/* End conditions */
c[1] = c[n] = 0.0;
/* Get cubic coefficients */
b[1] = (y[2] - y[1])/d[1] - d[i] * c[2];
c[1] = 0.0;
d[1] = c[2]/d[1];
b[n] = (y[n] - y[nm1])/d[nm1] + d[nm1] * c[nm1];
for(i=2 ; i<n ; i++) {
b[i] = (y[i+1]-y[i])/d[i] - d[i]*(c[i+1]+2.0*c[i]);
d[i] = (c[i+1]-c[i])/d[i];
c[i] = 3.0*c[i];
}
c[n] = 0.0;
d[n] = 0.0;
return;
}
/*
*	Splines a la Forsythe Malcolm and Moler
*	---------------------------------------
*	In this case the end-conditions are determined by fitting
*	cubic polynomials to the first and last 4 points and matching
*	the third derivitives of the spline at the end-points to the
*	third derivatives of these cubics at the end-points.
*/
static void
fmm_spline(int n, double *x, double *y, double *b, double *c, double *d)
{
int nm1, i;
double t;
/* Adjustment for 1-based arrays */
x--; y--; b--; c--; d--;
if(n < 2) {
errno = EDOM;
return;
}
if(n < 3) {
t = (y[2] - y[1]);
b[1] = t / (x[2]-x[1]);
b[2] = b[1];
c[1] = c[2] = d[1] = d[2] = 0.0;
return;
}
nm1 = n - 1;
/* Set up tridiagonal system */
/* b = diagonal, d = offdiagonal, c = right hand side */
d[1] = x[2] - x[1];
c[2] = (y[2] - y[1])/d[1];/* = +/- Inf	for x[1]=x[2] -- problem? */
for(i=2 ; i<n ; i++) {
d[i] = x[i+1] - x[i];
b[i] = 2.0 * (d[i-1] + d[i]);
c[i+1] = (y[i+1] - y[i])/d[i];
c[i] = c[i+1] - c[i];
}
/* End conditions. */
/* Third derivatives at x[0] and x[n-1] obtained */
/* from divided differences */
b[1] = -d[1];
b[n] = -d[nm1];
c[1] = c[n] = 0.0;
if(n > 3) {
c[1] = c[3]/(x[4]-x[2]) - c[2]/(x[3]-x[1]);
c[n] = c[nm1]/(x[n] - x[n-2]) - c[n-2]/(x[nm1]-x[n-3]);
c[1] = c[1]*d[1]*d[1]/(x[4]-x[1]);
c[n] = -c[n]*d[nm1]*d[nm1]/(x[n]-x[n-3]);
}
/* Gaussian elimination */
for(i=2 ; i<=n ; i++) {
t = d[i-1]/b[i-1];
b[i] = b[i] - t*d[i-1];
c[i] = c[i] - t*c[i-1];
}
/* Backward substitution */
c[n] = c[n]/b[n];
for(i=nm1 ; i>=1 ; i--)
c[i] = (c[i]-d[i]*c[i+1])/b[i];
/* c[i] is now the sigma[i-1] of the text */
/* Compute polynomial coefficients */
b[n] = (y[n] - y[n-1])/d[n-1] + d[n-1]*(c[n-1]+ 2.0*c[n]);
for(i=1 ; i<=nm1 ; i++) {
b[i] = (y[i+1]-y[i])/d[i] - d[i]*(c[i+1]+2.0*c[i]);
d[i] = (c[i+1]-c[i])/d[i];
c[i] = 3.0*c[i];
}
c[n] = 3.0*c[n];
d[n] = d[nm1];
return;
}
/*
*	Periodic Spline
*	---------------
*	The end conditions here match spline (and its derivatives)
*	at x[1] and x[n].
*
*	Note: There is an explicit check that the user has supplied
*	data with y[1] equal to y[n].
*/
static void
periodic_spline(int n, double *x, double *y, double *b, double *c, double *d)
{
double s;
int i, nm1;
// double *e = (double *) R_alloc(n, sizeof(double));
arma::vec ev(n);
double* e = ev.memptr();
/* Adjustment for 1-based arrays */
x--; y--; b--; c--; d--; e--;
if(n < 2 || y[1] != y[n]) {
errno = EDOM;
return;
}
if(n == 2) {
b[1] = b[2] = c[1] = c[2] = d[1] = d[2] = 0.0;
return;
} else if (n == 3) {
b[1] = b[2] = b[3] = -(y[1] - y[2])*(x[1] - 2*x[2] + x[3])/(x[3]-x[2])/(x[2]-x[1]);
c[1] = -3*(y[1]-y[2])/(x[3]-x[2])/(x[2]-x[1]);
c[2] = -c[1];
c[3] = c[1];
d[1] = -2*c[1]/3/(x[2]-x[1]);
d[2] = -d[1]*(x[2]-x[1])/(x[3]-x[2]);
d[3] = d[1];
return;
}
/* else --------- n >= 4 --------- */
nm1 = n-1;
/* Set up the matrix system */
/* A = diagonal	 B = off-diagonal  C = rhs */
#define A	b
#define B	d
#define C	c
B[1]  = x[2] - x[1];
B[nm1]= x[n] - x[nm1];
A[1] = 2.0 * (B[1] + B[nm1]);
C[1] = (y[2] - y[1])/B[1] - (y[n] - y[nm1])/B[nm1];
for(i = 2; i < n; i++) {
B[i] = x[i+1] - x[i];
A[i] = 2.0 * (B[i] + B[i-1]);
C[i] = (y[i+1] - y[i])/B[i] - (y[i] - y[i-1])/B[i-1];
}
/* Choleski decomposition */
#define L	b
#define M	d
#define E	e
L[1] = sqrt(A[1]);
E[1] = (x[n] - x[nm1])/L[1];
s = 0.0;
for(i = 1; i <= nm1 - 2; i++) {
M[i] = B[i]/L[i];
if(i != 1) E[i] = -E[i-1] * M[i-1] / L[i];
L[i+1] = sqrt(A[i+1]-M[i]*M[i]);
s = s + E[i] * E[i];
}
M[nm1-1] = (B[nm1-1] - E[nm1-2] * M[nm1-2])/L[nm1-1];
L[nm1] = sqrt(A[nm1] - M[nm1-1]*M[nm1-1] - s);
/* Forward Elimination */
#define Y	c
#define D	c
Y[1] = D[1]/L[1];
s = 0.0;
for(i=2 ; i<=nm1-1 ; i++) {
Y[i] = (D[i] - M[i-1]*Y[i-1])/L[i];
s = s + E[i-1] * Y[i-1];
}
Y[nm1] = (D[nm1] - M[nm1-1] * Y[nm1-1] - s) / L[nm1];
#define X	c
X[nm1] = Y[nm1]/L[nm1];
X[nm1-1] = (Y[nm1-1] - M[nm1-1] * X[nm1])/L[nm1-1];
for(i=nm1-2 ; i>=1 ; i--)
X[i] = (Y[i] - M[i] * X[i+1] - E[i] * X[nm1])/L[i];
/* Wrap around */
X[n] = X[1];
/* Compute polynomial coefficients */
for(i=1 ; i<=nm1 ; i++) {
s = x[i+1] - x[i];
b[i] = (y[i+1]-y[i])/s - s*(c[i+1]+2.0*c[i]);
d[i] = (c[i+1]-c[i])/s;
c[i] = 3.0*c[i];
}
b[n] = b[1];
c[n] = c[1];
d[n] = d[1];
return;
}
#undef A
#undef B
#undef C
#undef L
#undef M
#undef E
#undef Y
#undef D
#undef X
/* These were the public interfaces */
static void
spline_coef(int method, int n, double *x, double *y,
double *b, double *c, double *d)
{
switch(method) {
case 1:
periodic_spline(n, x, y, b, c, d);	break;
case 2:
natural_spline(n, x, y, b, c, d);	break;
case 3:
fmm_spline(n, x, y, b, c, d);	break;
}
}
static void
spline_eval(int method, int nu, double *u, double *v,
int n, double *x, double *y, double *b, double *c, double *d)
{
/* Evaluate  v[l] := spline(u[l], ...),	    l = 1,..,nu, i.e. 0:(nu-1)
* Nodes x[i], coef (y[i]; b[i],c[i],d[i]); i = 1,..,n , i.e. 0:(*n-1)
*/
const int n_1 = n - 1;
int i, j, k, l;
double ul, dx, tmp;
if(method == 1 && n > 1) { /* periodic */
dx = x[n_1] - x[0];
for(l = 0; l < nu; l++) {
v[l] = fmod(u[l]-x[0], dx);
if(v[l] < 0.0) v[l] += dx;
v[l] += x[0];
}
} else
for(l = 0; l < nu; l++) v[l] = u[l];
for(l = 0, i = 0; l < nu; l++) {
ul = v[l];
if(ul < x[i] || (i < n_1 && x[i+1] < ul)) {
/* reset i  such that  x[i] <= ul <= x[i+1] : */
i = 0;
j = n;
do {
k = (i+j)/2;
if(ul < x[k]) j = k;
else i = k;
} while(j > i+1);
}
dx = ul - x[i];
/* for natural splines extrapolate linearly left */
tmp = (method == 2 && ul < x[0]) ? 0.0 : d[i];
v[l] = y[i] + dx*(b[i] + dx*(c[i] + dx*tmp));
}
}
class SplineCoef {
public:
SplineCoef(arma::vec x, arma::vec y, int method = 2) : x(x), y(y), method(method) {
int n = x.n_elem;
// if(y.n_elem != x.n_elem) REprintf("inputs of different lengths");
b.resize(n); c.resize(n); d.resize(n);
b.zeros(); c.zeros(); d.zeros();
if (method==1 && y(0) != y(n-1))
y(n-1)=y(0);
spline_coef(method, x.n_elem, x.memptr(), y.memptr(), b.memptr(), c.memptr(), d.memptr());
}
arma::vec eval(arma::vec xout) {
int nu = xout.n_elem;
arma::vec yout(nu);
spline_eval(method, nu, xout.memptr(), yout.memptr(),
x.n_elem, x.memptr(), y.memptr(), b.memptr(), c.memptr(), d.memptr());
return yout;
}
double eval(double xout) {
double yout;
spline_eval(method, 1, &xout, &yout,
x.n_elem, x.memptr(), y.memptr(), b.memptr(), c.memptr(), d.memptr());
return yout;
}
static arma::vec eval(double xout, std::vector<SplineCoef>& z);
arma::vec x, y, b, c, d;
int method;
};
arma::vec SplineCoef::eval(double xout, std::vector<SplineCoef>& z)
{
arma::vec yout(z.size());
for (size_t i=0; i<z.size(); i++)
yout(i) = z[i].eval(xout);
return yout;
}

// namespace boost { namespace numeric { namespace odeint {
//       template <>
//       struct is_resizeable<arma::vec>
//       {
// 	typedef boost::true_type type;
// 	const static bool value = type::value;
//       };
//       template <>
//       struct same_size_impl<arma::vec, arma::vec>
//       {
// 	static bool same_size(const arma::vec& x, const arma::vec& y)
// 	{
// 	  return x.size() == y.size();
// 	}
//       };
//       template<>
//       struct resize_impl<arma::vec, arma::vec>
//       {
// 	static void resize(arma::vec& v1, const arma::vec& v2)
// 	{
// 	  v1.resize(v2.size());
// 	}
//       };
//     } } } // namespace boost::numeric::odeint
// typedef arma::vec state_type;
// typedef std::vector<state_type> vector_state_type;
// typedef std::vector<double> vector_times;
// struct push_back_state_and_time
// {
//     std::vector< state_type >& m_states;
//     std::vector< double >& m_times;
//     push_back_state_and_time( std::vector< state_type > &states , std::vector< double > &times )
//     : m_states( states ) , m_times( times ) { }
//     void operator()( const state_type &x , double t )
//     {
//         m_states.push_back( x );
//         m_times.push_back( t );
//     }
// };
// struct push_back_state
// {
//     std::vector< state_type >& m_states;
//     push_back_state( std::vector< state_type > &states)
//     : m_states( states ) { }
//     void operator()( const state_type &x , double t )
//     {
//         m_states.push_back( x );
//     }
// };
// struct Flow {
//   size_t from, to;
//   SplineCoef s;
//   bool use_log;
// };
// // old version of libarmadillo-dev: mat.as_col() not available
// arma::vec flatten(arma::mat m) {
//   arma::vec v(m.n_cols*m.n_rows);
//   for (size_t j=0; j<m.n_cols; j++)
//     for (size_t i=0; i<m.n_rows; i++)
//       v(j*m.n_rows+i) = m(i,j);
//   return v;
// }
// struct StateComponents {
//   arma::vec P, L;
// };
// struct StateComponentsCombined {
//   arma::vec times;
//   arma::mat P, L;
// };
// // forward declaration
// class Markov {
// public:
//   typedef std::vector<Flow> Flows;
//   Flows& flows;
//   double minTm;
//   StateComponents stateComponents(const state_type& x) {
//     arma::vec P = x(arma::span(0,nStates-1));
//     arma::vec L = x(arma::span(nStates,2*nStates-1));
//   }
//   Markov(Flows& flows, double minTm = 1.0e-8) : flows(flows), minTm(minTm) {
//     nGradients = nStates = 0;
//     for(size_t i=0; i<flows.size(); i++) {
//       nStates = std::max(std::max(nStates,flows[i].from),flows[i].to);
//     }
//     nStates++;
//   }
//   void operator() ( const state_type &x , state_type &dxdt , const double t )
//   {
//     arma::vec rates(flows.size());
//     arma::vec dPdt(nStates,arma::fill::zeros);
//     arma::vec dLdt;
//     // destructure the state vector
//     StateComponents c = stateComponents(x);
//     // P (transition probabilities)
//     for (size_t i=0; i<flows.size(); i++) {
//       rates(i) = flows[i].s.eval(flows[i].use_log ? log(t<minTm ? minTm : t) : t);
//       // if (flows[i].use_log) rates(i) = exp(rates(i));
//       double delta = c.P(flows[i].from)*rates(i);
//       dPdt(flows[i].to) += delta;
//       dPdt(flows[i].from) -= delta;
//     }
//     // L (length of stay)
//     dLdt = c.P;
//     for (size_t i=0, start=0; i<flows.size(); start+=flows[i].gradients.size(), i++) {
//       for (size_t j=0; j<flows[i].gradients.size(); j++) {
//       }
//     }
//     dxdt = arma::join_cols(dPdt,
// 			   arma::join_cols(dLdt,
//   }
// };
// StateComponentsCombined run(Markov &model, arma::vec p0, arma::vec times) {
//   using namespace boost::numeric::odeint;
//   vector_state_type states;
//   // times, outtimes and vtimes all have the same data
//   vector_times outtimes; // not strictly needed - we could use push_back_state()
//   vector_times vtimes = arma::conv_to< vector_times >::from(times);
//   state_type x = join_cols(p0, arma::zeros(model.nStates+2*model.nStates*model.nGradients));
//   // BOOST_STATIC_ASSERT( is_resizeable<state_type>::value == true );
//   integrate_times(make_dense_output( 1.0e-10 , 1.0e-10 , runge_kutta_dopri5< state_type >() ),
// 		  model, x,
// 		  vtimes.begin(),
// 		  vtimes.end(),
// 		  1.0,
// 		  push_back_state_and_time(states, outtimes));
//   size_t nx = states[0].size(), nTimes = times.size();
//   // combine the results
//   arma::mat combined(nTimes,nx);
//   for (size_t i=0; i<nTimes; i++)
//     combined(arma::span(i,i),arma::span(0,nx-1)) = states[i].t();
//   // destructure the combined matrix
//   arma::mat P = combined(arma::span::all,arma::span(0,model.nStates-1));
//   arma::mat L = combined(arma::span::all,arma::span(model.nStates,2*model.nStates-1));
// }
// // Rcpp::List runMarkovODE(arma::vec y0, arma::vec times, arma::vec tlam, arma::mat lam, Rcpp::List gradients, arma::ivec from, arma::ivec to, Rcpp::LogicalVector use_logs, double minTm = 1.0e-8) {
// RcppExport SEXP runMarkovODE(SEXP _y0, SEXP _times, SEXP _tlam, SEXP _lam, SEXP _gradients, SEXP _from, SEXP _to, SEXP _use_logs, SEXP _minTm) {
//   using namespace Rcpp;
//   arma::vec y0 = as<arma::vec>(_y0);
//   arma::vec times = as<arma::vec>(_times);
//   arma::vec tlam = as<arma::vec>(_tlam);
//   arma::mat lam = as<arma::mat>(_lam);
//   arma::ivec from = as<arma::ivec>(_from);
//   arma::ivec to = as<arma::ivec>(_to);
//   LogicalVector use_logs = as<LogicalVector>(_use_logs);
//   double minTm = as<double>(_minTm);
//   Markov::Flows flows;
//   for (size_t i=0; i<from.size(); i++) {
//     SplineCoef coef = SplineCoef(use_logs[i] ? log(tlam) : tlam, lam.col(i));
//     for (int j=0; j<gradients.size(); j++)
//                      bool(use_logs[i])});
//   }
//   Markov markov(flows, minTm);
//   auto report = run(markov, y0, times);
//   return wrap(List::create(Named("times")=report.times,
// 			   Named("P")=report.P,
// 			   Named("L")=report.L,
// }

RcppExport SEXP multistate_ddt(SEXP P_, SEXP Pu_, SEXP Q_, SEXP Qu_) {
using namespace arma;
// P(nstates,nobs), Pu(nstates,ncoef,nobs), Q(nstates,nstates,nobs), Qu(nstates*nstates,ncoef,nobs)
mat P = Rcpp::as<mat>(P_);
cube Pu = Rcpp::as<cube>(Pu_);
cube Q = Rcpp::as<cube>(Q_);
cube Qu = Rcpp::as<cube>(Qu_);
int nstates = P.n_rows;
int nobs=P.n_cols;
int ncoef=Pu.n_cols;
mat dPdt = P*0.0;   // nstates,nobs
cube dPudt = Pu*0.0; // nstates,ncoef,nobs
for(int i=0; i<nobs; i++) {
mat Qi = Q.slice(i);
rowvec Pi = conv_to<rowvec>::from(P.col(i));
vec rs = sum(Qi,1);
for (int j=0; j<nstates; j++)
Qi(j,j) = -rs(j);
dPdt.col(i) = conv_to<vec>::from(Pi*Qi); // is this valid?
cube Qui = Qu(span::all, span::all, span(i,i));
Qui.reshape(nstates,nstates,ncoef);
mat Pui = Pu.slice(i);
for (int j=0; j<ncoef; j++) {
vec rs = sum(Qui.slice(j),1);
for (int m=0; m<nstates; m++)
Qui(m,m,j) = -rs(m);
}
dPudt.slice(i) = (Pui.t() * Qi).t();
for(int j=0; j<ncoef; j++) {
rowvec PiQuij = Pi * Qui.slice(j);
for(int k=0; k<nstates; k++)
dPudt(k,j,i) += PiQuij(k);
}
}
return Rcpp::wrap(Rcpp::List::create(dPdt,dPudt));
}

// class ExpM {
// public:
//   arma::mat Qmat;
//   ExpM(arma::mat _Qmat) : Qmat(_Qmat) { }
//   void operator() ( const state_type &x , state_type &dxdt , const double t )
//   {
//     dxdt = (x.t() * Qmat).t();
//   }
// };
// RcppExport SEXP runExpM(SEXP _y0, SEXP _times, SEXP _Qmat) {
//   using namespace Rcpp;
//   arma::vec y0 = as<arma::vec>(_y0);
//   arma::vec times = as<arma::vec>(_times);
//   arma::mat Qmat = as<arma::mat>(_Qmat);
//   ExpM model(Qmat);
//   using namespace boost::numeric::odeint;
//   vector_state_type states;
//   // times, outtimes and vtimes all have the same data
//   vector_times outtimes; // not strictly needed - we could use push_back_state()
//   vector_times vtimes = arma::conv_to< vector_times >::from(times);
//   // BOOST_STATIC_ASSERT( is_resizeable<state_type>::value == true );
//   integrate_times(make_dense_output( 1.0e-10 , 1.0e-10 , runge_kutta_dopri5< state_type >() ),
// 		  model, y0,
// 		  vtimes.begin(),
// 		  vtimes.end(),
// 		  1.0,
// 		  push_back_state_and_time(states, outtimes));
//   size_t nx = states[0].size(), nTimes = times.size();
//   // combine the results
//   arma::mat combined(nTimes,nx);
//   for (size_t i=0; i<nTimes; i++)
//     combined(arma::span(i,i),arma::span(0,nx-1)) = states[i].t();
//   return wrap(combined);
// }
``````