Cran submission

.pre-commit-config.yaml

@@ -10,6 +10,7 @@ repos:
   hooks:
   - id: cppcheck
     args: ["--platform=win64", "--language=c++", "--std=c++17", "--enable=warning", "--quiet", "--inline-suppr", "--library=googletest", "--library=tinyxml2"]
+    files: \.cpp$
 
 - repo: https://github.com/pre-commit/mirrors-clang-format
   rev: v16.0.6

README.md

@@ -12,7 +12,7 @@ At its core, this package features multivariate spectral density estimation base
 
 Moreover, `dpca` implements the method to select the number of dynamic principal components of [Hallin & Liska (2007)](https://doi.org/10.1198/016214506000001275). Also the information criteria by [Bai & Ng, 2003](https://doi.org/10.1111/1468-0262.00273) are implemented, with the additional feature that the constant multplier in the penalty term is chosen in a data-driven way analogously to [Hallin & Liska (2007)](https://doi.org/10.1198/016214506000001275).
 
-## Things implemented so far:
+## Implementation roadmap:
 
 - [x] Estimation of multivariate spectral density using lag-window technique.
 - [x] Fast computation of dynamic eigenvalues/eigenvectors of spectrum using [`ARPACK`](https://en.wikipedia.org/wiki/ARPACK).
@@ -24,14 +24,15 @@ Moreover, `dpca` implements the method to select the number of dynamic principal
 - [ ] Forecasting methods.
 - [ ] Model assessment.
 - [ ] Port C code to modern C++.
+- [ ] Consider switch from ARPACK to [Spectra](https://spectralib.org/).
 
 We are aware of the R package [`freqdom`](https://CRAN.R-project.org/package=freqdom), developed by Siegfried Hörmann and Lukas Kidzinsiki which is a pure `R` implementation. `dpca` is written mainly in `C`. Although providing a similiar interface to that of `freqdom`, `dpca` has some unique features apart from being much faster.
 
 For instance, the convoluted filter which computes the dynamic common component from the output in `freqdom` is obtained by filtering the output twice: first to get the inputs \(what in `freqdom` is called "scores" in analogy to their FDA context\) and, second these inputs are filtered again to get the dynamic common component \("KLexpansion"\). `dpca` computes the convolution in the frequency domain. The advantage of this approach is that this filter is invariant with respect to multiplications of dynamic eigenvectors by a unit-length complex number.
 
 ## Installation
 
-`dpca` depends on [`ARPACK`](https://en.wikipedia.org/wiki/ARPACK) to compute dynamic eigenvalues/eigenvectors using the Implicitly Restarted Arnoldi Method which is much faster than R's `base::eigen()` whenever a truncated instead of the full spectral decomposition is required. `ARPACK` is shipped with `dpca`.
+`dpca` depends on [`ARPACK`](https://en.wikipedia.org/wiki/ARPACK) to compute dynamic eigenvalues/eigenvectors using the Implicitly Restarted Arnoldi Method which is much faster than R's `base::eigen()` whenever a truncated instead of the full spectral decomposition is required. ``dpca` ships `ARPACK`.
 
 To install `dpca` using `devtools`:
 

configure

@@ -1,25 +1,15 @@
-: ${R_HOME=`R RHOME`}
-if test -z "${R_HOME}"; then
-  echo "could not determine R_HOME"
-  exit 1
-fi
+# Copy required arpack sources into src/arpack
+# and configure with ISO C bindings
 
-CC=`"${R_HOME}/bin/R" CMD config CC`
-CFLAGS=`"${R_HOME}/bin/R" CMD config CFLAGS`
-CPPFLAGS=`"${R_HOME}/bin/R" CMD config CPPFLAGS`
-LDFLAGS=`"${R_HOME}/bin/R" CMD config LDFLAGS`
-FC=`"${R_HOME}/bin/R" CMD config FC`
-
-FSRC="z*.f icbazn.F90 debug.h stat.h"
-USRC="icnteq.f icopy.f iset.f iswap.f dvout.f ivout.f second_NONE.f zvout.f zmout.f"
+ZSRC="z*.f icbazn.F90 debug.h stat.h"
+USRC="icnteq.f icopy.f iset.f iswap.f second_NONE.f"
 FSRC_C="*.f icbazn.F90"
 
 tar -xf tools/arpack-ng-3.9.0.tar.gz -C tools/
-mkdir -p src/arpack
-for f in $FSRC; do cp tools/arpack-ng-3.9.0/SRC/$f src/arpack/; done
+
+for f in $ZSRC; do cp tools/arpack-ng-3.9.0/SRC/$f src/arpack/; done
+cat tools/arpack-ng-3.9.0/SRC/znaupd.f | sed -e '/if (msglvl .gt. 0) then/,/end if/s/^/c /' > src/arpack/znaupd.f
 for f in $USRC; do cp tools/arpack-ng-3.9.0/UTIL/$f src/arpack/; done
 cat tools/arpack-ng-3.9.0/arpackdef.h.in | sed -e 's/@INTERFACE64@/0/g' > src/arpack/arpackdef.h
 cat tools/arpack-ng-3.9.0/arpackicb.h.in | sed -e 's/@INTERFACE64@/0/g' > src/arpack/arpackicb.h
 cp tools/arpack-ng-3.9.0/ICB/arpack.h src/arpack/
-cd src/arpack &&
-  for f in $FSRC_C; do ${FC} -c -fpic $f; done

configure.win

@@ -1,26 +1,14 @@
-: ${R_HOME=`R RHOME`}
-if test -z "${R_HOME}"; then
-  echo "could not determine R_HOME"
-  exit 1
-fi
-
-CC=`"${R_HOME}/bin/R" CMD config CC`
-CFLAGS=`"${R_HOME}/bin/R" CMD config CFLAGS`
-CPPFLAGS=`"${R_HOME}/bin/R" CMD config CPPFLAGS`
-LDFLAGS=`"${R_HOME}/bin/R" CMD config LDFLAGS`
-FC=`"${R_HOME}/bin/R" CMD config FC`
+# Copy required arpack sources into src/arpack
+# and configure with ISO C bindings
 
 ZSRC="z*.f icbazn.F90 debug.h stat.h"
 USRC="icnteq.f icopy.f iset.f iswap.f dvout.f ivout.f second_NONE.f zvout.f zmout.f"
 FSRC_C="*.f icbazn.F90"
 
 tar -xf tools/arpack-ng-3.9.0.tar.gz -C tools/
-mkdir -p src/arpack
+
 for f in $ZSRC; do cp tools/arpack-ng-3.9.0/SRC/$f src/arpack/; done
 for f in $USRC; do cp tools/arpack-ng-3.9.0/UTIL/$f src/arpack/; done
 cat tools/arpack-ng-3.9.0/arpackdef.h.in | sed -e 's/@INTERFACE64@/0/g' > src/arpack/arpackdef.h
 cat tools/arpack-ng-3.9.0/arpackicb.h.in | sed -e 's/@INTERFACE64@/0/g' > src/arpack/arpackicb.h
 cp tools/arpack-ng-3.9.0/ICB/arpack.h src/arpack/
-cd src/arpack &&
-  for f in $FSRC_C; do ${FC} -c -fpic $f; done &&
-  ar cr libarpack.a *.o

src/Makevars

@@ -4,5 +4,14 @@
 # PKG_LIBS = -llapack -lblas -lgcov -L"arpack-ng-3.9.0/build" -l:libarpack.a
 # PKG_CFLAGS = -I"arpack-ng-3.9.0/ICB/" -I"arpack-ng-3.9.0/" -O0 --coverage
 # PKG_LIBS = -lblas -L"arpack-ng-3.9.0/build" -l:libarpack.a
-PKG_LIBS = arpack/*.o -lblas -llapack
+PKG_LIBS = -Larpack -larpack $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
 PKG_CFLAGS = -I"arpack" -DUSING_R
+
+$(SHLIB): arpack/libarpack.a
+
+arpack/libarpack.a:
+	@(cd arpack && $(MAKE) libarpack.a \
+		FC="$(FC)" FFLAGS="$(FFLAGS) $(FPICFLAGS)" AR="$(AR)")
+
+# PKG_LIBS = arpack/*.o $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+# PKG_CFLAGS = -I"arpack" -DUSING_R

src/Makevars.win

@@ -1,2 +1,8 @@
-PKG_LIBS = arpack/*.o  -llapack -lblas -lgfortran -lquadmath
+PKG_LIBS = -Larpack -larpack -llapack $(BLAS_LIBS) $(FLIBS)
 PKG_CFLAGS = -I"arpack" -DUSING_R
+
+$(SHLIB): arpack/libarpack.a
+
+arpack/libarpack.a:
+	@(cd arpack && $(MAKE) libarpack.a \
+		FC="$(FC)" FFLAGS="$(FFLAGS) $(FPICFLAGS)" AR="$(AR)")

src/arpack/Makefile

@@ -0,0 +1,18 @@
+ARPACK_SOURCES = $(wildcard *.f *.F90)
+ARPACK_OBJECTS = $(ARPACK_SOURCES:.f=.o)
+ARPACK_OBJECTS := $(ARPACK_OBJECTS:.F90=.o)
+ARPACK_FFLAGS = $(filter-out -pedantic, $(FFLAGS) $(FPICFLAGS))
+
+libarpack.a: $(ARPACK_OBJECTS)
+	$(AR) rcs $@ $^
+
+%.o: %.f
+	$(FC) $(ARPACK_FFLAGS) -c $< -o $@
+
+%.o: %.F90
+	$(FC) $(ARPACK_FFLAGS) -c $< -o $@
+
+clean:
+	rm -f *.o $(LIBRARY)
+
+.PHONY: clean

src/arpack/dvout.f

@@ -0,0 +1,6 @@
+      SUBROUTINE DVOUT(LOUT, N, SX, IDIGIT, IFMT)
+      CHARACTER*( *)    IFMT
+      INTEGER            IDIGIT, LOUT, N
+      DOUBLE PRECISION   SX( * )
+      RETURN
+      END

src/arpack/ivout.f

@@ -0,0 +1,5 @@
+      SUBROUTINE IVOUT(LOUT, N, IX, IDIGIT, IFMT)
+      INTEGER    IX(*), N, IDIGIT, LOUT
+      CHARACTER  IFMT*(*)
+      RETURN
+      END

src/arpack/zmout.f

@@ -0,0 +1,6 @@
+      SUBROUTINE ZMOUT(LOUT, M, N, A, LDA, IDIGIT, IFMT)
+      INTEGER            M, N, IDIGIT, LDA, LOUT
+      Complex*16         A( LDA, * )
+      CHARACTER          IFMT*( * )
+      RETURN
+      END

src/arpack/zvout.f

@@ -0,0 +1,6 @@
+      SUBROUTINE ZVOUT(LOUT, N, CX, IDIGIT, IFMT)
+      INTEGER            N, IDIGIT, LOUT
+      Complex*16         CX( * )
+      CHARACTER          IFMT*( * )
+      RETURN
+      END

src/dpca.c

@@ -78,7 +78,7 @@ SEXP R_dpca
     fourier_transform(covs, nrx, nrx, freqs, nfreqs, lags, nlags,
                       (double _Complex *)COMPLEX(spec));
 
-
+    SEXP hl_select;
     if (select_q) { // chose q ( number of dynamic factors) using Hallin & Liska (2007) method
 
         int lps = length(r_penalty_scales);
@@ -127,6 +127,21 @@ SEXP R_dpca
         R_Free(temp_evecs);
         R_Free(temp_evals);
 
+        hl_select = PROTECT(allocVector(VECSXP, 5));
+        SEXP chosen_q = PROTECT(allocVector(INTSXP, 1));
+        *INTEGER(chosen_q) = q;
+        SET_VECTOR_ELT(hl_select, 0, unpenalized_ic_vals);
+        SET_VECTOR_ELT(hl_select, 1, sample_var);
+        SET_VECTOR_ELT(hl_select, 2, chosen_q);
+        SET_VECTOR_ELT(hl_select, 3, info);
+        SET_VECTOR_ELT(hl_select, 4, q_path);
+        SEXP nms_hl_select = PROTECT(allocVector(STRSXP, 5));
+        SET_STRING_ELT(nms_hl_select, 0, mkChar("unpenalized_ic_vals"));
+        SET_STRING_ELT(nms_hl_select, 1, mkChar("sample_var"));
+        SET_STRING_ELT(nms_hl_select, 2, mkChar("q"));
+        SET_STRING_ELT(nms_hl_select, 3, mkChar("info"));
+        SET_STRING_ELT(nms_hl_select, 4, mkChar("q_path"));
+        setAttrib(hl_select, R_NamesSymbol, nms_hl_select);
     } else { // user-specfied q
         q = *INTEGER(r_q);
         evecs = PROTECT(alloc3DArray(CPLXSXP, q, nrx, nfreqs));
@@ -165,6 +180,7 @@ SEXP R_dpca
                 }
             }
         }
+        hl_select = PROTECT(allocVector(VECSXP, 0));
     }
 
     for (int i = 0; i < nfreqs; i++) {
@@ -211,27 +227,6 @@ SEXP R_dpca
     SET_STRING_ELT(nms_filter, 1, mkChar("dcc"));
     setAttrib(filter, R_NamesSymbol, nms_filter);
 
-    SEXP hl_select;
-    if (select_q) {
-        hl_select = PROTECT(allocVector(VECSXP, 5));
-        SEXP chosen_q = PROTECT(allocVector(INTSXP, 1));
-        *INTEGER(chosen_q) = q;
-        SET_VECTOR_ELT(hl_select, 0, unpenalized_ic_vals);
-        SET_VECTOR_ELT(hl_select, 1, sample_var);
-        SET_VECTOR_ELT(hl_select, 2, chosen_q);
-        SET_VECTOR_ELT(hl_select, 3, info);
-        SET_VECTOR_ELT(hl_select, 4, q_path);
-        SEXP nms_hl_select = PROTECT(allocVector(STRSXP, 5));
-        SET_STRING_ELT(nms_hl_select, 0, mkChar("unpenalized_ic_vals"));
-        SET_STRING_ELT(nms_hl_select, 1, mkChar("sample_var"));
-        SET_STRING_ELT(nms_hl_select, 2, mkChar("q"));
-        SET_STRING_ELT(nms_hl_select, 3, mkChar("info"));
-        SET_STRING_ELT(nms_hl_select, 4, mkChar("q_path"));
-        setAttrib(hl_select, R_NamesSymbol, nms_hl_select);
-    } else {
-        hl_select = PROTECT(allocVector(VECSXP, 0));
-    }
-
     SEXP res = PROTECT(allocVector(VECSXP, 7));
     SET_VECTOR_ELT(res, 0, spec);
     SET_VECTOR_ELT(res, 1, eig);

src/eigs.c

@@ -89,7 +89,6 @@ void arnoldi_eigs(double _Complex *mat, int dim, int ldm, int q,
 
   int cnt = 0;
   while (ido != 99) {
-    /* call arpack like you would have, but, use znaupd_c instead of znaupd_ */
     znaupd_c(&ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
              workd, workl, lworkl, rwork, &info);
 
@@ -106,7 +105,6 @@ void arnoldi_eigs(double _Complex *mat, int dim, int ldm, int q,
   }
 
 
-  /* call arpack like you would have, but, use zneupd_c instead of zneupd_ */
   zneupd_c(rvec, howmny, select, d, z, ldz, sigma, workev, bmat, N, which, nev,
            tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, rwork,
            &info);
@@ -164,3 +162,4 @@ void arnoldi_eigs(double _Complex *mat, int dim, int ldm, int q,
     }
   }
 }
+