From 36dd08c7c97689d2eb8c5fcacc48c46dc976528c Mon Sep 17 00:00:00 2001
From: daniel <daniel@kapla.at>
Date: Fri, 6 May 2022 22:28:08 +0200
Subject: [PATCH] forgot to add

---
 tensorPredictors/NAMESPACE  |  2 ++
 tensorPredictors/src/init.c |  6 +++++-
 tensorPredictors/src/ttm.c  | 18 +++++++++---------
 3 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/tensorPredictors/NAMESPACE b/tensorPredictors/NAMESPACE
index 5cff5c2..221aacf 100644
--- a/tensorPredictors/NAMESPACE
+++ b/tensorPredictors/NAMESPACE
@@ -13,6 +13,7 @@ export(approx.kronecker)
 export(colKronecker)
 export(dist.projection)
 export(dist.subspace)
+export(kpir.approx)
 export(kpir.base)
 export(kpir.kron)
 export(kpir.momentum)
@@ -20,6 +21,7 @@ export(kpir.new)
 export(mat)
 export(matpow)
 export(matrixImage)
+export(mcrossprod)
 export(reduce)
 export(rowKronecker)
 export(tensor_predictor)
diff --git a/tensorPredictors/src/init.c b/tensorPredictors/src/init.c
index 5640700..2a18f9a 100644
--- a/tensorPredictors/src/init.c
+++ b/tensorPredictors/src/init.c
@@ -9,10 +9,14 @@
 /* Tensor Times Matrix a.k.a. Mode Product */
 extern SEXP ttm(SEXP A, SEXP X, SEXP mode);
 
+/* Tensor Mode Covariance e.g. `(1 / n) * A_(m) A_(m)^T` */
+extern SEXP mcrossprod(SEXP A, SEXP mode);
+
 /* List of registered routines (e.g. C entry points) */
 static const R_CallMethodDef CallEntries[] = {
     // {"FastPOI_C_sub", (DL_FUNC) &FastPOI_C_sub, 5},  // NOT USED
-    {"C_ttm", (DL_FUNC) &ttm, 3},
+    {"C_ttm",        (DL_FUNC) &ttm,        3},
+    {"C_mcrossprod", (DL_FUNC) &mcrossprod, 2},
     {NULL, NULL, 0}
 };
 
diff --git a/tensorPredictors/src/ttm.c b/tensorPredictors/src/ttm.c
index 760dc43..51a24e1 100644
--- a/tensorPredictors/src/ttm.c
+++ b/tensorPredictors/src/ttm.c
@@ -11,7 +11,7 @@
 /**
  * Tensor Times Matrix a.k.a. Mode Product
  *
- * @param A multi-dimensionl array
+ * @param A multi-dimensional array
  * @param B matrix
  * @param m mode index (1-indexed)
  */
@@ -41,12 +41,12 @@ extern SEXP ttm(SEXP A, SEXP B, SEXP m) {
         error("Dimension missmatch (mode dim not equal to ncol)");
     }
 
-    // calc nr of response elements `prod(dim(X)[-mode]) * ncol(X)`,
-    int leny = 1;
+    // calc nr of response elements `prod(dim(A)[-mode]) * ncol(A)` (size of C),
+    int sizeC = 1;
     // and the strides
-    //  `stride[0] <- prod(dim(X)[seq_len(mode - 1)])`
-    //  `stride[1] <- dim(X)[mode]`
-    //  `stride[2] <- prod(dim(X)[-seq_len(mode)])`
+    //  `stride[0] <- prod(dim(A)[seq_len(mode - 1)])`
+    //  `stride[1] <- dim(A)[mode]`
+    //  `stride[2] <- prod(dim(A)[-seq_len(mode)])`
     int stride[3] = {1, INTEGER(dim)[mode], 1};
     for (int i = 0; i < length(dim); ++i) {
         int size = INTEGER(dim)[i];
@@ -54,7 +54,7 @@ extern SEXP ttm(SEXP A, SEXP B, SEXP m) {
         if (!size) {
             error("Zero dimension detected");
         }
-        leny *= (i == mode) ? nrows(B) : size;
+        sizeC *= (i == mode) ? nrows(B) : size;
         stride[0] *= (i < mode) ? size : 1;
         stride[2] *= (i > mode) ? size : 1;
     }
@@ -63,7 +63,7 @@ extern SEXP ttm(SEXP A, SEXP B, SEXP m) {
     int nrow = nrows(B);
 
     // create response object C
-    SEXP C = PROTECT(allocVector(REALSXP, leny));
+    SEXP C = PROTECT(allocVector(REALSXP, sizeC));
 
     // raw data access pointers
     double* a = REAL(A);
@@ -88,7 +88,7 @@ extern SEXP ttm(SEXP A, SEXP B, SEXP m) {
     }
     /*
     // Tensor Times Matrix / Mode Product (reference implementation)
-    memset(c, 0, leny * sizeof(double));
+    memset(c, 0, sizeC * sizeof(double));
     for (int i2 = 0; i2 < stride[2]; ++i2) {
         for (int i1 = 0; i1 < stride[1]; ++i1) { // stride[1] == ncols(B)
             for (int j = 0; j < nrow; ++j) {