Support different length scales for different dimensions in graph GPs…

… (ported from onnela-lab/gptools#23).

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: gptoolsStan
 Title: Gaussian Processes on Graphs and Lattices in 'Stan'
-Version: 0.1.0
+Version: 0.2.0
 Authors@R: c(
     person("Till", "Hoffmann", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-4403-0722")),

NEWS.md

@@ -1,3 +1,8 @@
 # gptoolsStan 0.1.0
 
 Initial release.
+
+# gptoolsStan 0.2.0
+
+- Graph-based approximations of Gaussian processes now support different length scales along different dimensions (ported from https://github.com/onnela-lab/gptools/pull/23). This change requires cmdstan 2.36.0 or above due to a bug in `gp_matern32_cov` (see https://github.com/stan-dev/math/pull/3084).
+- Functions ending `_log_abs_det_jacobian` have been renamed to `_log_abs_det_jac` to comply with upcoming changes (see https://github.com/stan-dev/stanc3/issues/1470 for details).

R/model.R

@@ -9,10 +9,26 @@
 #'
 #' # Compile the model with paths set up to include 'Stan' sources from 'gptoolsStan'.
 #' model <- cmdstan_model(
-#'   stan_file="/path/to/your/model.stan",
-#'   include_paths=gptools_include_path(),
+#'     stan_file = "/path/to/your/model.stan",
+#'     include_paths = gptools_include_path(),
 #' )
 #' }
 gptools_include_path <- function() {
-    system.file("extdata", package="gptoolsStan")
+    version <- package_version(cmdstanr::cmdstan_version())
+    if (version < package_version("2.36.0")) {
+        warning(
+            sprintf(
+                paste(
+                    "cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels",
+                    "with different length scales for different dimensions; see",
+                    "https://github.com/stan-dev/math/pull/3084 for details. Your",
+                    "model may yield unexpected results or crash if you use",
+                    "nearest-neighbor Gaussian processes with Matern 3/2 kernels.",
+                    "Your cmdstan version is %s; consider upgrading."
+                ),
+                version
+            )
+        )
+    }
+    system.file("extdata", package = "gptoolsStan")
 }

inst/extdata/gptools/fft1.stan

@@ -75,11 +75,14 @@ Evaluate the log absolute determinant of the Jacobian associated with
 :param n: Size of the real signal. Necessary because the size cannot be inferred from :code:`rfft`.
 :returns: Log absolute determinant of the Jacobian.
 */
-real gp_rfft_log_abs_det_jacobian(vector cov_rfft, int n) {
+real gp_rfft_log_abs_det_jac(vector cov_rfft, int n) {
     vector[n %/% 2 + 1] rfft_scale = gp_evaluate_rfft_scale(cov_rfft, n);
     return - sum(log(rfft_scale[1:n %/% 2 + 1])) -sum(log(rfft_scale[2:(n + 1) %/% 2]))
         - log(2) * ((n - 1) %/% 2) + n * log(n) / 2;
 }
+real gp_rfft_log_abs_det_jacobian(vector cov_rfft, int n) {
+    reject("`gp_rfft_log_abs_det_jacobian` has been renamed to `gp_rfft_log_abs_det_jac` to comply with upcoming changes (see https://github.com/stan-dev/stanc3/issues/1470 for details).");
+}
 
 
 /**
@@ -94,7 +97,7 @@ real gp_rfft_lpdf(vector y, vector loc, vector cov_rfft) {
     int n = size(y);
     int nrfft = n %/% 2 + 1;
     vector[n] z = gp_rfft(y, loc, cov_rfft);
-    return std_normal_lpdf(z) + gp_rfft_log_abs_det_jacobian(cov_rfft, n);
+    return std_normal_lpdf(z) + gp_rfft_log_abs_det_jac(cov_rfft, n);
 }
 
 

inst/extdata/gptools/fft2.stan

@@ -140,7 +140,7 @@ Evaluate the log absolute determinant of the Jacobian associated with
     Fourier transform :code:`cov_rfft2`).
 :returns: Log absolute determinant of the Jacobian.
 */
-real gp_rfft2_log_abs_det_jacobian(matrix cov_rfft2, int width) {
+real gp_rfft2_log_abs_det_jac(matrix cov_rfft2, int width) {
     int height = rows(cov_rfft2);
     int n = width * height;
     int fftwidth = width %/% 2 + 1;
@@ -171,6 +171,9 @@ real gp_rfft2_log_abs_det_jacobian(matrix cov_rfft2, int width) {
     ladj += - log2() * nterms + n * log(n) / 2;
     return ladj;
 }
+real gp_rfft2_log_abs_det_jacobian(matrix cov_rfft2, int width) {
+    reject("`gp_rfft2_log_abs_det_jacobian` has been renamed to `gp_rfft2_log_abs_det_jac` to comply with upcoming changes (see https://github.com/stan-dev/stanc3/issues/1470 for details).");
+}
 
 
 /**
@@ -185,7 +188,7 @@ Evaluate the log probability of a two-dimensional Gaussian process realization i
 */
 real gp_rfft2_lpdf(matrix y, matrix loc, matrix cov_rfft2) {
     return std_normal_lpdf(to_vector(gp_rfft2(y, loc, cov_rfft2)))
-        + gp_rfft2_log_abs_det_jacobian(cov_rfft2, cols(y));
+        + gp_rfft2_log_abs_det_jac(cov_rfft2, cols(y));
 }