Added Simon Wood's tensor product C backend (thanks! infinitely more

efficient), Zhenghua kindly modified the NOMAD naming to not conflict
with other conventions (thanks Brian Ripley for pointing this out).
commit. Version number bumped, uploaded to CRAN.

CHANGELOG

@@ -1,3 +1,11 @@
+Changes from Version 0.15-26 to 0.15-27 [06-May-2017]
+
+* Zhenguya fixed minor issue with dim.bs()
+
+* Incorporating Simon Wood's C-based tensor product code (was using his R-based code, the C code is substantially faster)
+
+* Naming of variables in NOMAD port aligned with Brian Ripley's comment ('Identifiers starting with an underscore followed by an upper-case letter or another underscore are reserved for system macros and should not be used in portable code (including not as guards in C/C++ headers).')
+
 Changes from Version 0.15-25 to 0.15-26 [29-April-2017]
 
 * Zhenghua patched up issues arising with windows build

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: crs
-Version: 0.15-26
-Date: 2017-05-02
+Version: 0.15-27
+Date: 2017-05-06
 Imports: boot, stats, np, quantreg, rgl
 Suggests: logspline, quadprog
 Title: Categorical Regression Splines

R/mgcv.R

@@ -1,21 +1,36 @@
+## Copyright (C) 2000-2012 Simon N. Wood  simon.wood@r-project.org
+## 
+## This program is free software; you can redistribute it and/or
+## modify it under the terms of the GNU General Public License
+## as published by the Free Software Foundation; either version 2
+## of the License, or (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## (www.gnu.org/copyleft/gpl.html)
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
+## USA.
 
-#############################################################
-## Type 1 tensor product methods start here (i.e. Wood, 2006)
-#############################################################
-
-tensor.prod.model.matrix<-function(X)
+tensor.prod.model.matrix <- function(X) {
 # X is a list of model matrices, from which a tensor product model matrix is to be produced.
 # e.g. ith row is basically X[[1]][i,]%x%X[[2]][i,]%x%X[[3]][i,], but this routine works 
-# column-wise, for efficiency
-{ m<-length(X)
-  X1<-X[[m]]
-  n<-nrow(X1)
-  if (m>1) for (i in (m-1):1)
-  { X0<-X1;X1<-matrix(0,n,0)
-    for (j in 1:ncol(X[[i]]))
-    X1<-cbind(X1,X[[i]][,j]*X0)
-  }
-  X1
+# column-wise, for efficiency, and does work in compiled code.
+  m <- length(X)              ## number to row tensor product
+  d <- unlist(lapply(X,ncol)) ## dimensions of each X
+  n <- nrow(X[[1]])           ## columns in each X
+  X <- as.numeric(unlist(X))  ## append X[[i]]s columnwise
+  T <- numeric(n*prod(d))     ## storage for result
+  .Call(mgcv_tmm,X,T,d,m,n)   ## produce product
+  ## Give T attributes of matrix. Note that initializing T as a matrix 
+  ## requires more time than forming the row tensor product itself (R 3.0.1)
+  attr(T,"dim") <- c(n,prod(d)) 
+  class(T) <- "matrix"
+  T
 } ## end tensor.prod.model.matrix
 
 uniquecombs<-function(x) {

R/zzz.R

@@ -1,5 +1,5 @@
 .onAttach <- function (lib, pkg) {
-	packageStartupMessage("Categorical Regression Splines (version 0.15-26)\n[vignette(\"crs_faq\") provides answers to frequently asked questions]", domain = NULL,  appendLF = TRUE)
+	packageStartupMessage("Categorical Regression Splines (version 0.15-27)\n[vignette(\"crs_faq\") provides answers to frequently asked questions]", domain = NULL,  appendLF = TRUE)
 }
 
 .onLoad <- function (lib, pkg) {

crsver

@@ -1 +1 @@
-0.15-26
+0.15-27

src/Makevars

@@ -1,6 +1,6 @@
 NOMAD_SRC           = ./nomad_src
 SGTELIB_SRC         = ./sgtelib_src
-OBJECTS             = RuniqueCombs.o bspline.o gsl_bspline.o snomadr.o crs_init.o $(NOMAD_SRC)/nomad.o \
+OBJECTS             = RuniqueCombs.o bspline.o gsl_bspline.o snomadr.o crs_init.o mgcv.o $(NOMAD_SRC)/nomad.o \
 											$(NOMAD_SRC)/Barrier.o $(NOMAD_SRC)/Cache.o $(NOMAD_SRC)/Cache_File_Point.o $(NOMAD_SRC)/Cache_Point.o \
 											$(NOMAD_SRC)/Cache_Search.o $(NOMAD_SRC)/Clock.o $(NOMAD_SRC)/Direction.o $(NOMAD_SRC)/Directions.o $(NOMAD_SRC)/Display.o \
 											$(NOMAD_SRC)/Double.o $(NOMAD_SRC)/Eval_Point.o $(NOMAD_SRC)/Evaluator.o $(NOMAD_SRC)/Evaluator_Control.o \

src/crs_init.c

@@ -2,6 +2,7 @@
 #include <Rinternals.h>
 #include <stdlib.h> // for NULL
 #include <R_ext/Rdynload.h>
+#include "mgcv.h"
 
 /* FIXME: 
    Check these declarations against the C/Fortran source code.
@@ -25,6 +26,7 @@ static const R_CMethodDef CEntries[] = {
 };
 
 static const R_CallMethodDef CallEntries[] = {
+    { "mgcv_tmm",         (DL_FUNC) &mgcv_tmm,          5}, 
     {"smultinomadRSolve", (DL_FUNC) &smultinomadRSolve, 1},
     {"snomadRInfo",       (DL_FUNC) &snomadRInfo,       1},
     {"snomadRSolve",      (DL_FUNC) &snomadRSolve,      1},

src/mgcv.c

@@ -0,0 +1,72 @@
+/* Copyright (C) 2000-2012 Simon N. Wood  simon.wood@r-project.org
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+(www.gnu.org/copyleft/gpl.html)
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
+USA. */
+
+#include "mgcv.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <R.h>
+
+void mgcv_tensor_mm(double *X,double *T,int *d,int *m,int *n) {
+/* Code for efficient production of row tensor product model matrices.
+   X contains rows of matrices to be producted. Contains m matrices,
+   d[i] is number of columns in ith matrix (which are stored in ascending 
+   order). Each column has n rows. T is the target matrix, with n rows
+   and \prod_i d[i] columns.   
+*/
+  ptrdiff_t start, i,j,k, tp=1, xp=0,pd;
+  double *Xj,*Xj1,*Xk, *Tk,*p,*p1,*p2;
+  /*Rprintf("m = %d  n = %d  d = ",*m,*n);
+    for (i=0;i<*m;i++) Rprintf(" %d,",d[i]);*/
+  /* compute total columns in X, xp, and total columns in T, tp ... */
+  for (i=0;i<*m;i++) { xp += d[i];tp *= d[i];} 
+  Xk = X + (xp-d[*m-1]) * *n; /* start of last matrix in X */
+  Tk = T + (tp-d[*m-1]) * *n; /* start of last (filled) block in T */
+  /* initialize by putting final matrix in X into end block of T... */
+  p = Xk;p1 = Xk + *n * (ptrdiff_t) d[*m-1];p2 = Tk;
+  for (;p<p1;p++,p2++) *p2 = *p;
+  pd = d[*m-1]; /* cols in filled (end) block of T */
+  for (i = *m - 2;i>=0;i--) { /* work down through matrices stored in X */
+    Xk -= *n * (ptrdiff_t) d[i]; /* start of ith X matrix */
+    Xj = Xk;
+    start = tp - pd * d[i]; /* start column of target block in T */
+    p = T + start * *n; /* start location in T */
+    for (j=0;j<d[i];j++) { /* work through columns of ith X matrix */
+      p1 = Tk; /* start of T block to be multiplied by jth col of ith X matrix*/
+      Xj1 = Xj + *n; /* end of col j of ith X matrix */
+      /* following multiplies end T block by jth col of ith X matrix, storing 
+         result in appropriate T block... */
+      for (k=0;k<pd;k++) for (p2 = Xj;p2<Xj1;p2++,p1++,p++) *p = *p1 * *p2;
+      Xj += *n; /* start of col j+1 of ith X matrix */
+    }
+    pd *= d[i]; /* number of cols in filled in T block */
+    Tk = T + (tp - pd) * *n; /* start of filled in T block */
+  }  
+} /* mgcv_tensor_mm */
+
+void mgcv_tmm(SEXP x,SEXP t,SEXP D,SEXP M, SEXP N) {
+  /* wrapper for calling mgcv_tensor_mm using .Call */
+  double *X,*T;
+  int *d,*m,*n;
+  X = REAL(x);T=REAL(t);
+  d = INTEGER(D);
+  m = INTEGER(M);
+  n = INTEGER(N);
+  mgcv_tensor_mm(X,T,d,m,n);
+}
+

src/mgcv.h

@@ -0,0 +1,22 @@
+/* Copyright (C) 2000-2012 Simon N. Wood  simon.wood@r-project.org
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+(www.gnu.org/copyleft/gpl.html)
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
+USA. */
+
+#include <Rinternals.h>
+#include <Rconfig.h>
+
+void mgcv_tmm(SEXP x,SEXP t,SEXP D,SEXP M, SEXP N);