B-splines / QR decomposition

[Mark Clements]

I have an initial implementation of B-splines for TMB based on the splines package (see below). Would this be useful to add to TMB?

I would also like to implement natural splines for TMB using B-splines and a QR decomposition of the second derivatives at the boundary knots (as per the ns() function). However, I could not find a TMB implementation of the QR decomposition. Can anyone suggest how to go about the QR decomposition, please?

Kindly, Mark.


library(TMB)
src <- "#include <TMB.hpp>      // Links in the TMB libraries
template<class Type>
struct spl_struct {
    int order,            /* order of the spline */
    ordm1,            /* order - 1 (3 for cubic splines) */
    nknots,            /* number of knots */
    curs,            /* current position in knots vector */
    boundary;        /* must have knots[curs] <= x < knots[curs+1] */
                /* except for the boundary case */
  vector<Type> ldel;          /* differences from knots on the left */
  vector<Type> rdel;            /* differences from knots on the right */
  vector<Type> knots;            /* knot vector */
  vector<Type> coeff;            /* coefficients */
  vector<Type> a;            /* scratch array */
};
template<class Type>
int
set_cursor(spl_struct<Type> *sp, Type x)
{
    int i;
    /* don't assume x's are sorted */
    sp->curs = -1; /* Wall */
    sp->boundary = 0;
    for (i = 0; i < sp->nknots; i++) {
    if (sp->knots[i] >= x) sp->curs = i;
    if (sp->knots[i] > x) break;
    }
    if (sp->curs > sp->nknots - sp->order) {
    int lastLegit = sp->nknots - sp->order;
    if (x == sp->knots[lastLegit]) {
        sp->boundary = 1; sp->curs = lastLegit;
    }
    }
    return sp->curs;
}
template<class Type>
void
diff_table(spl_struct<Type> *sp, Type x, int ndiff)
{
  int i;
  for (i = 0; i < ndiff; i++) {
      sp->rdel[i] = sp->knots[sp->curs + i] - x;
      sp->ldel[i] = x - sp->knots[sp->curs - (i + 1)];
  }
}
/* fast evaluation of basis functions */
template<class Type>
void
basis_funcs(spl_struct<Type> *sp, Type x, matrix<Type> & b, int offset)
{
  diff_table(sp, x, sp->ordm1);
  b(0,offset) = Type(1);
  for (int j = 1; j <= sp->ordm1; j++) {
    Type saved = Type(0);
    for (int r = 0; r < j; r++) { // do not divide by zero
      Type den = sp->rdel[r] + sp->ldel[j - 1 - r];
      if(den != Type(0)) {
    Type term = b(r,offset)/den;
    b(r,offset) = saved + sp->rdel[r] * term;
    saved = sp->ldel[j - 1 - r] * term;
      } else {
    if(r != 0 || sp->rdel[r] != Type(0))
      b(r,offset) = saved;
    saved = Type(0);
      }
    }
    b(j,offset) = saved;
  }
}
template<class Type>
Type
evaluate(spl_struct<Type> *sp, Type x, int nder)
{
  // Type *lpt, *rpt, *apt, *ti = sp->knots + sp->curs;
  int ti = sp->curs,
    lpt, apt, rpt, inner,
    outer = sp->ordm1;
  if (sp->boundary && nder == sp->ordm1) { /* value is arbitrary */
    return Type(0);
  }
  while(nder--) {  // FIXME: divides by zero
    // for(inner = outer, apt = sp->a, lpt = ti - outer; inner--; apt++, lpt++)
    for(inner = outer, apt = 0, lpt = ti - outer; inner--; apt++, lpt++)
      //*apt = outer * (*(apt + 1) - *apt)/(*(lpt + outer) - *lpt);
      sp->a[apt] = Type(outer) * (sp->a[apt + 1] - sp->a[apt])/(sp->knots[lpt + outer] - sp->knots[lpt]);
    outer--;
  }
  diff_table(sp, x, outer);
  while(outer--)
    //for(apt = sp->a, lpt = sp->ldel + outer, rpt = sp->rdel, inner = outer + 1;
    for(apt = 0, lpt = outer, rpt = 0, inner = outer + 1;
        inner--; lpt--, rpt++, apt++)
        // FIXME: divides by zero
        sp->a[apt] = (sp->a[apt + 1] * sp->ldel[lpt] + sp->a[apt] * sp->rdel[rpt])/(sp->rdel[rpt] + sp->ldel[lpt]);
    return sp->a[0];
}
/* called from    splineDesign() : */
template<class Type>
matrix<Type>
spline_basis(vector<Type> kk,
         int ord,
         vector<Type> xx,
         vector<int> ders)
{
/* evaluate the non-zero B-spline basis functions (or their derivatives)
 * at xvals.  */
  int nk = kk.size();
  int nx = xx.size();
  int nd = ders.size();
  spl_struct<Type> * sp = new spl_struct<Type>();
  /* fill sp : */
  sp->order = ord;
  sp->ordm1 = ord - 1;
  sp->rdel = vector<Type>(sp->ordm1);
  sp->ldel = vector<Type>(sp->ordm1);
  sp->knots = kk; sp->nknots = nk;
  sp->a = vector<Type>(sp->order);
  matrix<Type> valM(sp->order, nx);
  vector<int> ioff(nx);
  for(int i = 0; i < nx; i++) {
    set_cursor(sp, xx[i]);
    int io = ioff[i] = sp->curs - sp->order;
    if (io < 0 || io > nk) {
      for (int j = 0; j < sp->order; j++) {
    valM(j,i) = Type(0); // R_NaN;
      }
    } else if (ders[i % nd] > 0) { /* slow method for derivatives */
      for(int ii = 0; ii < sp->order; ii++) {
    for(int j = 0; j < sp->order; j++) sp->a[j] = Type(0);
    sp->a[ii] = Type(1);
    valM(ii,i) =
      evaluate(sp, xx[i], ders[i % nd]);
      }
    } else {        /* fast method for value */
      basis_funcs(sp, xx[i], valM, i);
    }
  }
  delete sp;
  return valM.transpose();
}
template<class Type>
Type objective_function<Type>::operator() ()
{
  PARAMETER_VECTOR(beta);
  DATA_VECTOR(kk);
  DATA_INTEGER(ord);
  DATA_IVECTOR(ders);
  DATA_VECTOR(xx);
  matrix<Type> X = spline_basis(kk, ord, xx, ders);
  REPORT(X);
  return sum(X*beta);
}
"
write(src, file="splines.cpp")
compile("splines.cpp")
p <- list(beta=c(0,1,2,3))
dyn.load(dynlib("splines"))
f = MakeADFun(data=list(xx=c(1.1,1.2),kk=c(1,1,1,1,2,2,2,2),ord=4L, ders=0L), parameters=p,DLL="splines")
f$report(f$par)$X
## comparison
require(splines)
bs(c(1.1,1.2),Boundary.knots=c(1,2),df=4,intercept=TRUE)

[Hans Skaug]

There is some spline functionality in the file ./TMB/inst/include/tmbutils/splines.hpp on github,

but apparently not for B-splines. When I have used B-splines myself in the past,

I have set them up in R, and read in the coefficient matrix into TMB/ADMB and estimated

spline coefficients as parameters.

Regarding QR decomposition, it may be possible to access it through Eigen:

https://eigen.tuxfamily.org/dox/group__QR__Module.html

in the same way as you can access eigen decomposition:

https://github.com/kaskr/adcomp/issues/144

[Kasper Kristensen]

Hi Mark,

I'd suggest you share your code using the TMB:::install.contrib, see

https://github.com/kaskr/adcomp/wiki/User-contributed-code

On Tuesday, June 27, 2017 at 5:53:42 PM UTC+2, Mark Clements wrote:

[Jeff Laake]

I agree with Hans. What is useful about that approach is that lets you use any of the spline capabilities with other R formula terms.  Essentially you pass the design matrix to TMB containing the spline and other terms for the predictor variables. I've done this in my marked package on CRAN (also github) for TMB and ADMB and also for mark-recapture analysis in MARK with my RMark package.  This approach also works with implementing random effects via formula. It simplifies end-user software because they only need to know the R formula notation and any additional components for splines.

--jeff

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+unsubscribe@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/a74eac7f-b88d-4cc5-adaa-6ebbcd583c47%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Hans Skaug]

Yes, and you can easily do penalized splines:

1. Set up a B-spline basis in R and import the design matrix into TMB

2. Let spline coefficients be random effects. Their variance is the penalty parameter that gets estimated by maximum likelihood

3. You can let the spline coeffs follow 1st or 2nd order random walks to impose more smoothness.

When you do a ADREPORT() on the resulting spline function, the uncertainty in the penalty (variance) feeds into the uncertainty of the spline,

which I think it should.

[Mark Clements]

We have several examples where we need splines at the C++ level.

First, the motivating example is for accelerated failure time models, where we model survival at time t given covariates x as S(t|x)=S_0(t exp(x'beta)). Note that the evaluation points for S_0() will change with revised betas. We want to use a flexible parametric model for S_0() using a smooth function s_0() (e.g. natural splines), where S(t|x)=exp(-exp(s_0(log(t)+x'beta))). As a consequence, the design matrix for s_0() is not constant.

Note that the proportional hazards model S(t|x)=exp(-exp(s_0(log(t))+x'beta)) would have a constant design matrix for s_0() - see our package rstpm2 on CRAN.

Second, we have used splines to define log-transition hazards for a smooth Markov model using ordinary differential equations (e.g. https://github.com/mclements/purged/).

Note that the B-splines are useful for natural splines, P-splines and M-splines/I-splines.

I now realise that the QR matrix for the natural splines is being evaluated at the boundary knots and is hence constant, so I can calculate that Q matrix in R and pass it to TMB. I will try to package this up as per the suggested mechanism.

Finally: great work on TMB. This is an amazingly useful package.

Kindly, Mark.

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/1d08799c-fcdf-4989-8c23-e26f42a37ac5%40googlegroups.com.

[Mark Clements]

In re-implementing the splines::bs() function in TMB, I have run into a
weird bug (reproducible example attached). When I use the objective
function motivated by an accelerated failure time model:

template<class Type>
Type objective_function<Type>::operator() ()
{

PARAMETER(beta);
PARAMETER_VECTOR(betas);
DATA_VECTOR(x0);
DATA_VECTOR(x1);
DATA_VECTOR(boundaryKnots);
DATA_VECTOR(interiorKnots);
bs<Type> s1(boundaryKnots,interiorKnots,1);
vector<Type> xs = x0 + x1 * beta;
matrix<Type> Xs = s1.basis(xs);
REPORT(Xs);
Type f = (Xs*betas).sum();
return f;
}

for some given inputs, I get the correct spline basis (Xs) and function
value at the initial values. When I evaluate at a different value, I
also get the correct spline basis (ok), however I get an incorrect
objective function. This is particularly strange, as the objective is
only the sum of the basis matrix (which is correct) times the betas
vector...

Re-coding this with RcppEigen gives the correct objective function with
no evidence of a memory leak from valgrind (code available).

Is my implementation compatible with TMB? Any help with this would be
appreciated - I have been struggling with this for a while.

Kindly, Mark.

[Mark Clements]

To possibly answer my own question:
https://github.com/kaskr/adcomp/wiki/Things-you-should-NOT-do-in-TMB

The spline knots are external and do not depend on the parameters, while
the spline evaluation points and spline basis do depend on the
parameters. In summary, splines with accelerated failure time models
probably do not fit within the current TMB framework. Out of interest,
would this work in ADMB?

As an aside, extracting the Q matrix from a QR decomposition is possible
with TMB. The main trick was to #include <R_ext/Applic.h> before
#include <TMB.hpp>. The code is:

MatrixXd qr_q(const MatrixXd& X, double tol = 1E-12)
{
// Initialize member data and allocate heap memory
int n=X.rows(), p=X.cols(), rank=0;
MatrixXd qr = X, y(n,n), q(n,n);
int* pivot=(int*)R_alloc(p,sizeof(int));
double* tau=(double*)R_alloc(p,sizeof(double));
double* work=(double*)R_alloc(p*2,sizeof(double));
for(int i=0;i<p;i++)
pivot[i]=i+1;
y.setIdentity();
// LINPACK QR factorization via householder transformations
F77_CALL(dqrdc2)(qr.data(), &n, &n, &p, &tol, &rank, tau, pivot,
work);
// Compute orthogonal factor Q
F77_CALL(dqrqy)(qr.data(), &n, &rank, tau, y.data(), &n, q.data());
return q;
}

Kindly, Mark.

[Kasper Kristensen]

You are right that the problem is caused by parameter dependent branching. ADMB does not have this issue because it re-builds the computational graph for each function evaluation. You can rebuild the computational graph in TMB as well using 'obj$retape()' (after setting obj$env$parameters) within each evaluation, but it will get very slow so is not recommended.

Other options to deal with parameter dependent branching in TMB:

1. Use conditional expressions. E.g. CppAD::CondExpGt(x, y, expr_true, expr_false) returns 'expr_true' if x>y and 'expr_false' otherwise. variants 'Gt' (>), 'Ge' (>=), 'Lt' (<), 'Le' (<=), 'Eq' (==) exist. However, this approach would require major re-write of the spline code. Especially all parts depending on the 'cursor'.

2. In principle one could implement an 'atomic function' that automatically re-builds the computational graph of the spline for each function evaluation (which is much cheaper than rebuilding the entire graph). This is in principle possible but there is not yet an example that demonstrates it.

________________________________________
From: tmb-...@googlegroups.com [tmb-...@googlegroups.com] on behalf of Mark Clements [mark.c...@ki.se]
Sent: Tuesday, July 25, 2017 8:42 PM
To: Hans Skaug; TMB Users
Subject: Re: [TMB users] Re: B-splines / QR decomposition

Kindly, Mark.

--

To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/6E62DCB74E59D246B410609BB7A9D723018F24EC53%40KIMSX02.user.ki.se.