Merge pull request #7 from nlmixrdevelopment/issue-5

Issue 5

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: lotri
 Title: A Simple Way to Specify Symmetric, Block Diagonal Matrices
-Version: 0.2.2
+Version: 0.3.1
 Authors@R: c(person("Matthew L.","Fidler",
              role = c("aut", "cre"), email = "[email protected]",
              comment=c(ORCID="0000-0001-8538-6691")))

NEWS.md

@@ -1,5 +1,7 @@
-# lotri XXX
+# lotri 0.3.1
 * Change errors/warnings to use `call.=FALSE` or equivalent. 
+* Refactor C code to reduce complexity
+* Change C code to play nicely with `rchk`
 # lotri 0.2.2
 * Bug fix for conditional matrices
 * Now accessing `$lower` and `$upper` gives default values even if it

README.Rmd

@@ -146,7 +146,7 @@ To do this, you simply have to enclose the properties after the
 conditional variable.  That is `et1 ~ id(lower=3)`.
 
 
-## Combining symmetric named matrices
+## Combining symmetric (named) matrices
 
 Now there is even a faster way to do a similar banded matrix
 concatenation with `lotriMat`
@@ -182,6 +182,33 @@ print(mb)
 autoplot(mb)
 ```
 
+You may also combine named and unnamed matrices, but the resulting
+matrix will be unnamed, and still be faster than `Matrix`:
+
+```{r}
+testList <- list(lotri({et2 + et3 + et4 ~ c(40,
+                            0.1, 20,
+                            0.1, 0.1, 30)}),
+                     lotri(et5 ~ 6),
+                     lotri(et1+et6 ~c(0.1, 0.01, 1)),
+                     matrix(c(1L, 0L, 0L, 1L), 2, 2))
+
+matf <- function(.mats){
+  .omega <- as.matrix(Matrix::bdiag(.mats))
+  return(.omega)
+}
+
+print(matf(testList))
+
+print(lotriMat(testList))
+
+mb <- microbenchmark(matf(testList),lotriMat(testList))
+
+print(mb)
+
+autoplot(mb)
+```
+
 ## New features
 
 A new feature is the ability to condition on variables by `|`.  This

README.md

@@ -192,7 +192,7 @@ To do this, you simply have to enclose the properties after the
 conditional variable.  That is `et1 ~ id(lower=3)`.
 
 
-## Combining symmetric named matrices
+## Combining symmetric (named) matrices
 
 Now there is even a faster way to do a similar banded matrix
 concatenation with `lotriMat`
@@ -244,16 +244,69 @@ mb <- microbenchmark(matf(testList),lotriMat(testList))
 
 print(mb)
 #> Unit: microseconds
-#>                expr     min       lq      mean   median       uq      max neval
-#>      matf(testList) 464.486 470.1440 554.46817 474.5055 519.6730 6347.862   100
-#>  lotriMat(testList)   1.109   1.4095   2.59312   1.9420   2.5075   59.198   100
+#>                expr     min      lq      mean  median       uq      max neval
+#>      matf(testList) 463.599 491.021 598.57232 521.231 585.1135 5710.686   100
+#>  lotriMat(testList)   1.297   1.646   2.56874   2.732   3.0105    7.734   100
 
 autoplot(mb)
 #> Coordinate system already present. Adding new coordinate system, which will replace the existing one.
 ```
 
 <img src="man/figures/README-unnamed-chunk-8-1.png" title="plot of chunk unnamed-chunk-8" alt="plot of chunk unnamed-chunk-8" width="100%" />
 
+You may also combine named and unnamed matrices, but the resulting
+matrix will be unnamed, and still be faster than `Matrix`:
+
+
+```r
+testList <- list(lotri({et2 + et3 + et4 ~ c(40,
+                            0.1, 20,
+                            0.1, 0.1, 30)}),
+                     lotri(et5 ~ 6),
+                     lotri(et1+et6 ~c(0.1, 0.01, 1)),
+                     matrix(c(1L, 0L, 0L, 1L), 2, 2))
+
+matf <- function(.mats){
+  .omega <- as.matrix(Matrix::bdiag(.mats))
+  return(.omega)
+}
+
+print(matf(testList))
+#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
+#> [1,] 40.0  0.1  0.1    0 0.00 0.00    0    0
+#> [2,]  0.1 20.0  0.1    0 0.00 0.00    0    0
+#> [3,]  0.1  0.1 30.0    0 0.00 0.00    0    0
+#> [4,]  0.0  0.0  0.0    6 0.00 0.00    0    0
+#> [5,]  0.0  0.0  0.0    0 0.10 0.01    0    0
+#> [6,]  0.0  0.0  0.0    0 0.01 1.00    0    0
+#> [7,]  0.0  0.0  0.0    0 0.00 0.00    1    0
+#> [8,]  0.0  0.0  0.0    0 0.00 0.00    0    1
+
+print(lotriMat(testList))
+#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
+#> [1,] 40.0  0.1  0.1    0 0.00 0.00    0    0
+#> [2,]  0.1 20.0  0.1    0 0.00 0.00    0    0
+#> [3,]  0.1  0.1 30.0    0 0.00 0.00    0    0
+#> [4,]  0.0  0.0  0.0    6 0.00 0.00    0    0
+#> [5,]  0.0  0.0  0.0    0 0.10 0.01    0    0
+#> [6,]  0.0  0.0  0.0    0 0.01 1.00    0    0
+#> [7,]  0.0  0.0  0.0    0 0.00 0.00    1    0
+#> [8,]  0.0  0.0  0.0    0 0.00 0.00    0    1
+
+mb <- microbenchmark(matf(testList),lotriMat(testList))
+
+print(mb)
+#> Unit: nanoseconds
+#>                expr    min       lq      mean   median       uq     max neval
+#>      matf(testList) 448951 454581.5 535816.43 463802.5 528361.0 2228159   100
+#>  lotriMat(testList)    989   1281.5   2690.28   2396.0   2731.5   26770   100
+
+autoplot(mb)
+#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.
+```
+
+<img src="man/figures/README-unnamed-chunk-9-1.png" title="plot of chunk unnamed-chunk-9" alt="plot of chunk unnamed-chunk-9" width="100%" />
+
 ## New features
 
 A new feature is the ability to condition on variables by `|`.  This

src/asLotriMat.c

@@ -20,7 +20,7 @@ SEXP _asLotriMatGen(SEXP x, SEXP extra, SEXP def, SEXP dims, SEXP dimn, const ch
     UNPROTECT(pro);
     Rf_errorcall(R_NilValue, "extra properties need default try 'lotri(matrix,x=3,default=\"id\")'");
   }
-  SEXP extraNames = Rf_getAttrib(extra, R_NamesSymbol);
+  SEXP extraNames = PROTECT(Rf_getAttrib(extra, R_NamesSymbol)); pro++;
   for (int i = lExtra; i--;) {
     if (Rf_isNull(VECTOR_ELT(extra, i))) {
       nNull++;

src/lotriBounds.c

@@ -48,9 +48,9 @@ SEXP _lotriGetBounds(SEXP lst_, SEXP format, SEXP startNum) {
     Rf_errorcall(R_NilValue, _("when format is specified, 'startNum' must be a single integer"));
   }
   int len = Rf_length(li.lst);
-  int totdim=0;
+  int totdim=0, named = 1;
   for (int i = 0; i < len; ++i) {
-    totdim += getCheckDim(li.lst, i);
+    totdim += getCheckDim(li.lst, i, &named);
   }
   SEXP retN = PROTECT(Rf_allocVector(STRSXP, totdim)); pro++;
   SEXP boundLower = PROTECT(Rf_allocVector(REALSXP, totdim)); pro++;

src/lotriLstToMat.c

@@ -27,10 +27,10 @@ SEXP lotriToLstMat(SEXP lotri) {
   return ret;
 }
 
-lotriInfo assertCorrectMatrixProperties(SEXP lst_, SEXP format, SEXP startNum) {
+lotriInfo assertCorrectMatrixProperties(SEXP lst_, SEXP format, SEXP startNum, int *named) {
   int type = TYPEOF(lst_);
   if (type != VECSXP) {
-    if (isSymNameMat(lst_)) {
+    if (isSymNameMat(lst_, *named)) {
       lotriInfo li;
       li.sym = 1;
       li.lst = R_NilValue;
@@ -54,17 +54,18 @@ lotriInfo assertCorrectMatrixProperties(SEXP lst_, SEXP format, SEXP startNum) {
 
 SEXP _lotriLstToMat(SEXP lst_, SEXP format, SEXP startNum) {
   int type, totN, pro = 0;
-  lotriInfo li = assertCorrectMatrixProperties(lst_, format, startNum);
+  int named = 2;
+  lotriInfo li = assertCorrectMatrixProperties(lst_, format, startNum, &named);
   if (li.sym) return lst_;
   PROTECT(li.lst); pro++;
   int len = Rf_length(li.lst);
   int totdim = 0;
-  int i, j;
+  int i;
   if (len == 2) {
     int repN = isSingleInt(VECTOR_ELT(li.lst, 1), NA_INTEGER);
     if (repN == NA_INTEGER){
     } else if (repN > 0) {
-      if (isSymNameMat(VECTOR_ELT(li.lst, 0))){
+      if (isSymNameMat(VECTOR_ELT(li.lst, 0), named)){
 	SEXP new = PROTECT(Rf_allocVector(VECSXP, 1)); pro++;
 	SET_VECTOR_ELT(new, 0, li.lst);
 	SEXP ret = _lotriLstToMat(new, format, startNum);
@@ -74,7 +75,7 @@ SEXP _lotriLstToMat(SEXP lst_, SEXP format, SEXP startNum) {
     }
   }
   for (i = 0; i < len; ++i) {
-    totdim += getCheckDim(li.lst, i);
+    totdim += getCheckDim(li.lst, i, &named);
   }
   SEXP ret = PROTECT(Rf_allocMatrix(REALSXP, totdim, totdim)); pro++;
   SEXP retN = PROTECT(Rf_allocVector(STRSXP, totdim)); pro++;
@@ -83,8 +84,6 @@ SEXP _lotriLstToMat(SEXP lst_, SEXP format, SEXP startNum) {
   memset(retd, 0, sizeof(double)*totdim*totdim);
   // Now use memcpy/ integer conversion to c
   SEXP cur;
-  double *curd;
-  int *curi;
   int curBand = 0;
   SEXP dimnames, colnames, sameS;
   int nsame;
@@ -99,38 +98,19 @@ SEXP _lotriLstToMat(SEXP lst_, SEXP format, SEXP startNum) {
       type = TYPEOF(cur);
     }
     totN = Rf_ncols(cur);
-    dimnames = Rf_getAttrib(cur, R_DimNamesSymbol);
-    colnames = VECTOR_ELT(dimnames, 1);
-    for (int cursame = nsame; cursame--;){
-      if (type == REALSXP) {
-	curd = REAL(cur);
-	for (j = 0; j  < totN; ++j) {
-	  memcpy(&retd[totdim*(curBand+j)+curBand],
-		 &curd[totN*j], sizeof(double)*totN);
-	  // Repeats dim names of repeated matrices
-	  setStrElt(retN, colnames, curBand, j,
-		    li.fmt, li.doFormat, &li.counter, nsame);
-	}
-      } else {
-	curi = INTEGER(cur);
-	for (j = 0; j < totN; ++j) {
-	  double *to = &retd[totdim*(curBand+j)+curBand];
-	  double *last = to + totN; // N - count
-	  int *from = &curi[totN*j];
-	  while (to != last) {
-	    *(to++) = (double)(*(from++));
-	  }
-	  setStrElt(retN, colnames, curBand, j,
-		    li.fmt, li.doFormat, &li.counter, nsame);
-	}
-      }
-      curBand += totN;
+    if (named) {
+      dimnames = Rf_getAttrib(cur, R_DimNamesSymbol);
+      colnames = VECTOR_ELT(dimnames, 1);
     }
+    lotriLstToMatFillInMatrixBand(retd, nsame, type, named, totN, totdim,
+				  retN, colnames, &curBand, &li, cur);
+  }
+  if (named) {
+    dimnames = PROTECT(Rf_allocVector(VECSXP, 2)); pro++;
+    SET_VECTOR_ELT(dimnames, 0, retN);
+    SET_VECTOR_ELT(dimnames, 1, retN);
+    Rf_setAttrib(ret, R_DimNamesSymbol, dimnames);
   }
-  dimnames = PROTECT(Rf_allocVector(VECSXP, 2)); pro++;
-  SET_VECTOR_ELT(dimnames, 0, retN);
-  SET_VECTOR_ELT(dimnames, 1, retN);
-  Rf_setAttrib(ret, R_DimNamesSymbol, dimnames);
   UNPROTECT(pro);
   return ret;
 }

src/lotriLstToMat.h

@@ -13,7 +13,7 @@ typedef struct lotriInfo {
 
 SEXP lotriToLstMat(SEXP lotri);
 
-lotriInfo assertCorrectMatrixProperties(SEXP lst_, SEXP format, SEXP startNum);
+lotriInfo assertCorrectMatrixProperties(SEXP lst_, SEXP format, SEXP startNum, int *named);
 SEXP _lotriLstToMat(SEXP lst_, SEXP format, SEXP startNum);
 
 static inline lotriInfo _lotriLstToMat0(SEXP lst_, SEXP format, SEXP startNum) {
@@ -55,4 +55,38 @@ static inline lotriInfo _lotriLstToMat0(SEXP lst_, SEXP format, SEXP startNum) {
   return ret;
 }
 
+static inline void lotriLstToMatFillInMatrixBand(double *retd, int nsame, int type, int named, int totN, int totdim,
+						 SEXP retN, SEXP colnames, int *curBand, lotriInfo *li,
+						 SEXP cur) {
+  for (int cursame = nsame; cursame--;){
+    if (type == REALSXP) {
+      double *curd = REAL(cur);
+      for (int j = 0; j  < totN; ++j) {
+	memcpy(&retd[totdim*(*curBand+j)+(*curBand)],
+	       &curd[totN*j], sizeof(double)*totN);
+	// Repeats dim names of repeated matrices
+	if (named) {
+	  setStrElt(retN, colnames, (*curBand), j,
+		    li->fmt, li->doFormat, &(li->counter), nsame);
+	}
+      }
+    } else {
+      int *curi = INTEGER(cur);
+      for (int j = 0; j < totN; ++j) {
+	double *to = &retd[totdim*(*curBand+j)+(*curBand)];
+	double *last = to + totN; // N - count
+	int *from = &curi[totN*j];
+	while (to != last) {
+	  *(to++) = (double)(*(from++));
+	}
+	if (named) {
+	  setStrElt(retN, colnames, (*curBand), j,
+		    li->fmt, li->doFormat, &(li->counter), nsame);
+	}
+      }
+    }
+    *curBand += totN;
+  }
+}
+
 #endif

src/matlist.h

@@ -40,13 +40,14 @@ static inline int isSingleInt(SEXP in, int defaultVal) {
   return defaultVal;
 }
 
-static inline int isSymNameMat(SEXP cur) {
+static inline int isSymNameMat(SEXP cur, int named) {
   int type = TYPEOF(cur);
   if (type == INTSXP || type == REALSXP) {
     if (Rf_isMatrix(cur)){
       int nrows = Rf_nrows(cur);
       int ncols = Rf_ncols(cur);
       if (nrows == ncols) {
+	if (!named) return nrows;
 	SEXP dimn = Rf_getAttrib(cur, R_DimNamesSymbol);
 	if (dimn != R_NilValue) {
 	  return nrows;
@@ -57,7 +58,7 @@ static inline int isSymNameMat(SEXP cur) {
   return 0;
 }
 
-static inline int getCheckDim(SEXP lst, int i) {
+static inline int getCheckDim(SEXP lst, int i, int *named) {
   SEXP cur = VECTOR_ELT(lst, i);
   int type = TYPEOF(cur);
   int same=1;
@@ -75,11 +76,20 @@ static inline int getCheckDim(SEXP lst, int i) {
     cur = VECTOR_ELT(cur, 0);
     type = TYPEOF(cur);
   }
-  int ret = isSymNameMat(cur);
+  int ret = isSymNameMat(cur, *named);
   if (ret){
     return ret*same;
   } else {
-    Rf_errorcall(R_NilValue, _("list element %d is not a symmetric named matrix"), i+1);
+    // if named is 2, then reassign named to 0 and return the dimension, reset the named to 0
+    if (*named == 2) {
+      ret = isSymNameMat(cur, 0);
+      if (ret) {
+	*named = 0;
+	return ret*same;
+      }
+    }
+    if (*named) Rf_errorcall(R_NilValue, _("list element %d is not a symmetric named matrix"), i+1);
+    else Rf_errorcall(R_NilValue, _("list element %d is not a symmetric matrix"), i+1);
   }
   return 0;
 }