https://github.com/cran/rstpm2
Tip revision: 936fa1723eb80b0f9d4b12a1e50e5c542b549b35 authored by Mark Clements on 03 March 2021, 16:10:02 UTC
version 1.5.2
version 1.5.2
Tip revision: 936fa17
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));
}