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 > × )
// : 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;
// std::vector<SplineCoef> gradients;
// 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;
// arma::mat gradP, gradL;
// };
// struct StateComponentsCombined {
// arma::vec times;
// arma::mat P, L;
// arma::cube gradP, gradL;
// };
// // forward declaration
// class Markov {
// public:
// typedef std::vector<Flow> Flows;
// Flows& flows;
// double minTm;
// size_t nStates, nGradients;
// 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));
// arma::mat gradP(arma::vec(x(arma::span(2*nStates,2*nStates+nStates*nGradients-1))).memptr(),
// nStates, nGradients);
// arma::mat gradL(arma::vec(x(arma::span(2*nStates+nStates*nGradients,
// 2*nStates+2*nStates*nGradients-1))).memptr(),
// nStates, nGradients);
// return {P, L, gradP, gradL};
// }
// Markov(Flows& flows, double minTm = 1.0e-8) : flows(flows), minTm(minTm) {
// nGradients = nStates = 0;
// for(size_t i=0; i<flows.size(); i++) {
// nGradients += flows[i].gradients.size();
// 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;
// arma::mat dGradPdt(nStates,nGradients,arma::fill::zeros);
// arma::mat dGradLdt;
// // 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;
// // gradients for 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++) {
// double delta = c.gradP(flows[i].from,start+j)*rates(i)+
// c.P(flows[i].from)*flows[i].gradients[j].eval(t);
// dGradPdt(flows[i].to,start+j) += delta;
// dGradPdt(flows[i].from,start+j) -= delta;
// }
// }
// // gradients for L
// dGradLdt = c.gradP;
// dxdt = arma::join_cols(dPdt,
// arma::join_cols(dLdt,
// arma::join_cols(flatten(dGradPdt),
// flatten(dGradLdt))));
// }
// };
// 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));
// arma::cube gradP(arma::mat(combined(arma::span::all,
// arma::span(2*model.nStates,2*model.nStates+model.nStates*model.nGradients-1))).memptr(),
// nTimes, model.nStates, model.nGradients);
// arma::cube gradL(arma::mat(combined(arma::span::all,
// arma::span(2*model.nStates+model.nStates*model.nGradients,
// 2*model.nStates+2*model.nStates*model.nGradients-1))).memptr(),
// nTimes, model.nStates, model.nGradients);
// return {times, P, L, gradP, gradL};
// }
// // 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);
// List gradients = as<List>(_gradients);
// 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));
// std::vector<SplineCoef> vgradients;
// arma::mat g = Rcpp::as<arma::mat>(gradients[i]);
// for (int j=0; j<gradients.size(); j++)
// vgradients.push_back(SplineCoef(tlam, g.col(j)));
// flows.push_back({size_t(from(i)), size_t(to(i)), coef, vgradients,
// 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,
// Named("gradP")=report.gradP,
// Named("gradL")=report.gradL));
// }
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);
// }
