fixed coupla_defination.stan of zero_constrain

examples/linear_algebra/R/correlation_angles_example.R

@@ -1,12 +1,15 @@
 library(cmdstanr)
 fp <- file.path("./examples/linear_algebra/stan/correlation_angles_example.stan")
-mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra")
+mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra", force_recompile = T)
 
 T <- matrix(c(0.8, -0.9, -0.9,
               -0.9, 1.1, 0.3,
               -0.9, 0.4, 0.9), 3, 3)
 
-
+T <- matrix(c(1, -0.9, -0.9,
+              -0.9, 1, 0.3,
+              -0.9, 0.4, 1), 3, 3)
+chol(T)
 mod_out <- mod$optimize(
   data = list(N = 3,
               R = T)
@@ -17,12 +20,10 @@ mod_out <- mod$sample(
               R = T),
   chains = 2,
   seed = 23421,
+
   parallel_chains = 2,
   iter_warmup = 600,
   iter_sampling = 600
 )
 
-mod_out$summary("R_hat")
-mat <- matrix(mod_out$summary("R_hat")$mean, 3, 3)
-mat
-chol(mat)
+mod_out$summary()

examples/linear_algebra/R/correlation_angles_example2.R

@@ -0,0 +1,120 @@
+library(cmdstanr)
+fp <- file.path("./examples/linear_algebra/stan/correlation_angles_constrain_example.stan")
+mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra", force_recompile = T)
+
+library(rstan)
+expose_stan_functions(fp)
+
+library(StanHeaders)
+stanFunction("multiply_lower_tri_self_transpose", x = L)
+
+
+N <- 3
+test_mat <- matrix(c(0.8, -0.9, -0.9,
+              -0.9, 1.1, 0.3,
+              -0.9, 0.4, 0.9), 3, 3)
+
+T <- matrix(c(1, -0.9, -0.9,
+              -0.9, 1, 0.3,
+              -0.9, 0.4, 1), 3, 3)
+
+
+test_mat <- matrix(c(1.0000,   0,         0,         0,   -0.9360,
+                     0,    1.0000,   -0.5500,   -0.3645,   -0.5300,
+                     0,   -0.5500,    1.0000,   0,    0.0875,
+                     0,   -0.3645,    0,    1.0000,    0.4557,
+                     -0.9360,   -0.5300,    0.0875,    0.4557,    1.0000),
+                   byrow = T, 5, 5)
+
+where_zero <- (test_mat[upper.tri(test_mat)] == 0) * 1
+sum(where_zero) * 1
+
+fit <- stan(file = fp, 
+            data = list(N = nrow(test_mat),
+                        R = test_mat,
+                        is_symmetric = 1,
+                        Z = sum(where_zero),
+                        K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+                        where_zero = where_zero ),
+            chains = 1)
+
+K <- 5
+theta_test <- build_angle_mat(where_zero, runif(6, 0, pi), K)
+L = zero_constrain(theta_test, K)
+tcrossprod(L)
+
+angle_raw <-  runif(6)
+N <- length(where_zero);
+mat <- matrix(0, N, N)
+count <- 1;
+raw_count <- 1;
+
+for (k in 2:K){
+  mat[k, k] = 0;
+  for (j in 1:k - 1) {
+    if (where_zero[raw_count] != 1) {
+      mat[k, j] = angle_raw[count];
+      count <= count +  1;
+    }
+    raw_count = raw_count + 1;
+  }
+}
+
+count <- 0
+
+for (k in 2:K){
+  for (j in 1:(k - 1)) {
+    count <- count + 1
+    print(count)
+  }}
+
+
+
+mod_out <- mod$sample(
+  data = list(N = nrow(test_mat),
+              R = test_mat,
+              is_symmetric = 1,
+              Z = sum(where_zero),
+              K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+              where_zero = where_zero ),
+  chains = 2,
+  #init = 0,
+  #seed = 23421,
+  #adapt_delta = 0.8,
+  #max_treedepth = 10,
+  parallel_chains = 4,
+  iter_warmup = 400,
+  iter_sampling = 400
+)
+
+N <- nrow(test_mat)
+
+round(matrix(mod_out$summary("R_out")$mean,nrow(test_mat), nrow(test_mat)), 3)
+
+
+chol(test_mat)
+
+mod_out <- mod$optimize(
+  data = list(N = nrow(test_mat),
+                      R = test_mat,
+                      is_symmetric = 1,
+                      Z = sum(where_zero),
+                      K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+                      where_zero = where_zero ))
+
+round(matrix(mod_out$summary("R_out")$estimate,nrow(test_mat), nrow(test_mat)), 3)
+
+mod_out <- mod$sample(
+  data = list(N = nrow(test_mat),
+              R = test_mat,
+              is_symmetric = 1),
+  chains = 2,
+  seed = 23421,
+  adapt_delta = 0.9,
+  max_treedepth = 10,
+  parallel_chains = 2,
+  iter_warmup = 200,
+  iter_sampling = 200
+)
+mod_out$summary("R_out")
+round(matrix(mod_out$summary("R_out")$mean, N, N), 3)

examples/linear_algebra/stan/correlation_angles_constrain.stan

@@ -0,0 +1,100 @@
+functions {
+#include correlation_angles.stan
+#include triangular.stan
+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        out[counter] = M[i, j];
+        counter += 1;
+      }
+    }
+    return out;
+  }
+
+   vector lower_elements_constrain(matrix M, int A){
+    int n = rows(M);
+    int counter = 1;
+    vector[A] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        if(M[i, j] > 0){
+         out[counter] = M[i, j];
+         counter += 1;
+        }
+      }
+    }
+    return out;
+  }
+
+ matrix build_angle_mat (vector where_zero, vector angle_raw, int K) {
+  int N = num_elements(where_zero);
+  matrix[K, K] mat;
+  int count = 1;
+  int raw_count = 0;
+
+  mat[1, 1] = 0.0;
+  for (k in 2:K){
+    mat[k, k] = 0.0;
+    for (j in 1:k - 1) {
+      raw_count += 1;
+      if (where_zero[raw_count] != 1) {
+         mat[k, j] = angle_raw[count];
+         count += 1;
+      }
+    }
+  }
+
+  return mat;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    if (i > 2) {
+      for (j in 2:(i - 1)) {
+        real prod_sines = prod(sin(angle[i, 1:j - 1]));
+        real cos_theta;
+        if (is_nan(angle_raw[i, j]) == 1) {
+          cos_theta = -(dot_product(inv_mat[j, 1:j - 1], inv_mat[i, 1:j - 1]) ) / ( inv_mat[j, j] * prod_sines );
+            if ( cos_theta < -1 || cos_theta > 1 ) reject("cos_theta is ", cos_theta, " and must be in [-1, 1]"); // cos_theta = 0;
+           // else if( cos_theta > 1) cos_theta = 1;
+            angle[i, j] = acos( cos_theta );
+          }
+        inv_mat[i, j] = cos(angle[i, j]) * prod_sines;
+      }
+    }
+    inv_mat[i, i] = prod(sin(angle[i, 1:(i - 1)]));
+  }
+  return inv_mat;
+}

examples/linear_algebra/stan/correlation_angles_constrain.stanfunctions

@@ -0,0 +1,100 @@
+functions {
+#include correlation_angles.stan
+#include triangular.stan
+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        out[counter] = M[i, j];
+        counter += 1;
+      }
+    }
+    return out;
+  }
+
+   vector lower_elements_constrain(matrix M, int A){
+    int n = rows(M);
+    int counter = 1;
+    vector[A] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        if(M[i, j] > 0){
+         out[counter] = M[i, j];
+         counter += 1;
+        }
+      }
+    }
+    return out;
+  }
+
+ matrix build_angle_mat (vector where_zero, vector angle_raw, int K) {
+  int N = num_elements(where_zero);
+  matrix[K, K] mat;
+  int count = 1;
+  int raw_count = 0;
+
+  mat[1, 1] = 0.0;
+  for (k in 2:K){
+    mat[k, k] = 0.0;
+    for (j in 1:k - 1) {
+      raw_count += 1;
+      if (where_zero[raw_count] != 1) {
+         mat[k, j] = angle_raw[count];
+         count += 1;
+      }
+    }
+  }
+
+  return mat;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    if (i > 2) {
+      for (j in 2:(i - 1)) {
+        real prod_sines = prod(sin(angle[i, 1:j - 1]));
+        real cos_theta;
+        if (is_nan(angle_raw[i, j]) == 1) {
+          cos_theta = -(dot_product(inv_mat[j, 1:j - 1], inv_mat[i, 1:j - 1]) ) / ( inv_mat[j, j] * prod_sines );
+            if ( cos_theta < -1 || cos_theta > 1 ) reject("cos_theta is ", cos_theta, " and must be in [-1, 1]"); // cos_theta = 0;
+           // else if( cos_theta > 1) cos_theta = 1;
+            angle[i, j] = acos( cos_theta );
+          }
+        inv_mat[i, j] = cos(angle[i, j]) * prod_sines;
+      }
+    }
+    inv_mat[i, i] = prod(sin(angle[i, 1:(i - 1)]));
+  }
+  return inv_mat;
+}