diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6a48fb8..ef774aa 100644 --- a/.pre-commit-config.yaml +++ b/.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 diff --git a/README.md b/README.md index ed3a81c..0d33c17 100644 --- a/README.md +++ b/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,6 +24,7 @@ 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. @@ -31,7 +32,7 @@ For instance, the convoluted filter which computes the dynamic common component ## 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`: diff --git a/configure b/configure index e81306d..f2f35cc 100755 --- a/configure +++ b/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 diff --git a/configure.win b/configure.win index 116abb5..4905f89 100755 --- a/configure.win +++ b/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 diff --git a/src/Makevars b/src/Makevars index 0d8fe60..25cecae 100644 --- a/src/Makevars +++ b/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 diff --git a/src/Makevars.win b/src/Makevars.win index 8ef7572..0322bb7 100644 --- a/src/Makevars.win +++ b/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)") diff --git a/src/arpack/Makefile b/src/arpack/Makefile new file mode 100644 index 0000000..b5c4e7f --- /dev/null +++ b/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 diff --git a/src/arpack/dvout.f b/src/arpack/dvout.f new file mode 100644 index 0000000..98a86a4 --- /dev/null +++ b/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 diff --git a/src/arpack/ivout.f b/src/arpack/ivout.f new file mode 100644 index 0000000..ccc7872 --- /dev/null +++ b/src/arpack/ivout.f @@ -0,0 +1,5 @@ + SUBROUTINE IVOUT(LOUT, N, IX, IDIGIT, IFMT) + INTEGER IX(*), N, IDIGIT, LOUT + CHARACTER IFMT*(*) + RETURN + END diff --git a/src/arpack/zmout.f b/src/arpack/zmout.f new file mode 100644 index 0000000..f960af4 --- /dev/null +++ b/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 diff --git a/src/arpack/zvout.f b/src/arpack/zvout.f new file mode 100644 index 0000000..ebd135f --- /dev/null +++ b/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 diff --git a/src/dpca.c b/src/dpca.c index f90b331..6ce3559 100644 --- a/src/dpca.c +++ b/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); diff --git a/src/eigs.c b/src/eigs.c index 1d22aed..ccbb59b 100644 --- a/src/eigs.c +++ b/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, } } } +