diff --git a/tmp-save/maybe-fortuto-qcss.R b/tmp-save/maybe-fortuto-qcss.R
new file mode 100644
index 0000000..bda5472
--- /dev/null
+++ b/tmp-save/maybe-fortuto-qcss.R
@@ -0,0 +1,68 @@
+library(testthat)
+library(bigsnpr)
+
+bedfile <- file.path(tempdir(), "tmp-data/public-data3.bed")
+if (!file.exists(rdsfile <- sub_bed(bedfile, ".rds"))) {
+  zip <- tempfile(fileext = ".zip")
+  download.file(
+    "https://github.com/privefl/bigsnpr/blob/master/data-raw/public-data3.zip?raw=true",
+    destfile = zip, mode = "wb")
+  unzip(zip, exdir = tempdir())
+  rds <- snp_readBed(bedfile)
+  expect_identical(normalizePath(rds), normalizePath(rdsfile))
+}
+
+obj.bigSNP <- snp_attach(rdsfile)
+G <- obj.bigSNP$genotypes
+y <- obj.bigSNP$fam$affection
+POS2 <- obj.bigSNP$map$genetic.dist + 1000 * obj.bigSNP$map$chromosome
+
+sumstats <- bigreadr::fread2(file.path(tempdir(), "tmp-data/public-data3-sumstats.txt"))
+sumstats$n_eff <- sumstats$N
+map <- setNames(obj.bigSNP$map[-3], c("chr", "rsid", "pos", "a1", "a0"))
+df_beta <- snp_match(sumstats, map, join_by_pos = FALSE)
+
+ind_var <- df_beta$`_NUM_ID_`
+corr0 <- snp_cor_extendedThr(
+  G, ind.col = ind_var, thr_r2 = 0.2,
+  infos.pos = POS2[ind_var], size = 1 / 1000, ncores = 2)
+corr <- as_SFBM(corr0)
+
+maf <- snp_MAF(G, ind.col = ind_var)
+sd0 <- sqrt(2 * maf * (1 - maf))
+
+expect_equal(
+  snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, max_run = 2),
+  snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, max_run = 2)
+)
+#
+# expect_equal(
+#   snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, max_run = 1),
+#   snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, max_run = 1, ncores = 2)
+# )
+# 1 iter is fine
+# 2 iter is fine with ncores = 1
+# any is fine with ncores = 1
+# the problem comes from ncores = 2 and max_run > 1
+
+# test0 <- replicate(10, snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, max_run = 2))
+# #18378 / 10129 / 19061 is TRUE
+
+expect_equal(
+  test1 <- snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, max_run = 2, ncores = 3),
+  test2 <- snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, max_run = 2, ncores = 3)
+)
+
+# only r2_max is the same
+# thr are all 0.2 / keep are all TRUE
+
+expect_equal(
+  test3 <- snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, ncores = 3),
+  test4 <- snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, ncores = 3)
+)
+
+plot(lengths(test1$ind_imp), lengths(test2$ind_imp))
+
+res <- snp_qcimp_sumstats(corr, df_beta, sd0 = sd0, ncores = 2)
+expect_false(any(res$rm_qc, na.rm = TRUE))
+expect_equal(snp_qcimp_sumstats(corr, df_beta,  sd0 = sd0, ncores = 2), res)
diff --git a/tmp-tests/test-gwas-combi.R b/tmp-tests/test-gwas-combi.R
new file mode 100644
index 0000000..0621e35
--- /dev/null
+++ b/tmp-tests/test-gwas-combi.R
@@ -0,0 +1,21 @@
+library(bigsnpr)
+obj <- snp_attachExtdata()
+G <- obj$genotypes
+set.seed(1)
+y <- rnorm(nrow(G))
+w <- runif(5); w <- w / sum(w)
+G_w <- G[, 1:5] %*% w
+
+sapply(1:5, function(j) {
+  summary(lm(y ~ G[, j]))$coef[2, 1:3]
+})
+
+gwas <- big_univLinReg(G, y, ind.col = 1:5)
+gwas$score
+
+(Z_w <- crossprod(gwas$score, w))
+summary(mylm <- lm(y ~ G_w))$coef  # 0.570 instead of 0.266
+# crossprod(scale(G_w), scale(y)) / (sd(mylm$residuals) * sqrt(nrow(G)))
+
+var <- big_colstats(G, ind.col = 1:5)$var
+(beta_w <- crossprod(gwas$estim * var, w) / var(G_w))  # both 0.0739
diff --git a/tmp-tests/test-gwas-multipop.R b/tmp-tests/test-gwas-multipop.R
new file mode 100644
index 0000000..d8233a5
--- /dev/null
+++ b/tmp-tests/test-gwas-multipop.R
@@ -0,0 +1,137 @@
+library(bigsnpr)
+
+celiac <- snp_attach("../Dubois2010_data/FinnuncorrNLITUK1UK3hap300_QC_norel.rds")
+dim(G <- celiac$genotypes)  # 15155 x 281122
+CHR <- celiac$map$chromosome
+POS <- celiac$map$physical.pos
+y <- celiac$fam$affection - 1
+NCORES <- nb_cores()
+
+pop <- snp_getSampleInfos(celiac, "../Dubois2010_data/FinnNLITUK1UK3.clusterv2")[[1]]
+pop2 <- c("Netherlands", "Italy", "UK1", "UK2", "Finland")[pop]
+table(pop2, exclude = NULL)
+# Finland       Italy Netherlands         UK1         UK2
+#    2436        1035        1623        3325        6736
+
+big_counts(G, ind.col = 1:10)
+G2 <- snp_fastImputeSimple(G, method = "mean2", ncores = NCORES)
+
+obj.svd <- snp_autoSVD(G2, CHR, POS, ncores = NCORES)
+plot(obj.svd)
+plot(obj.svd, type = "scores", scores = 1:6, coef = 0.6)
+PC <- predict(obj.svd)[, 1:6]
+
+pop_centers <- bigutilsr::geometric_median(PC, by_grp = pop2)
+
+ldist_to_centers <- apply(pop_centers, 1, function(center) {
+  log(rowSums(sweep(PC, 2, center, '-')^2))
+})
+hist(ldist_to_centers[pop2 != "Finland", "Finland"])
+hist(ldist_to_centers[pop2 == "Italy", "Italy"])
+hist(ldist_to_centers[pop2 != "Italy", "Italy"])
+hist(ldist_to_centers[pop2 != "UK2", "UK2"])
+
+ind_NW <- which(ldist_to_centers[, "UK2"] < 6)
+ind_Fin <- which(ldist_to_centers[, "Finland"] < 7)
+ind_It <- which(ldist_to_centers[, "Italy"] < 6)
+
+run_GWAS <- function(ind) {
+  big_univLinReg(G2, y.train = y[ind], ind.train = ind,
+                 covar.train = cbind(celiac$fam$sex, PC)[ind, ],
+                 ncores = NCORES)
+}
+
+gwas_NW <- run_GWAS(ind_NW)
+gwas_Fin <- run_GWAS(ind_Fin)
+gwas_It <- run_GWAS(ind_It)
+
+plot_grid(
+  snp_manhattan(gwas_NW,  CHR, POS, npoints = 20e3) +
+    ggplot2::scale_y_log10(),
+  snp_manhattan(gwas_Fin, CHR, POS, npoints = 20e3) +
+    ggplot2::scale_y_log10(),
+  snp_manhattan(gwas_It,  CHR, POS, npoints = 20e3) +
+    ggplot2::scale_y_log10(),
+  ncol = 1
+)
+
+ind_sig <- which(predict(gwas_NW, log10 = FALSE) < 1e-20)
+plot(gwas_NW$estim[ind_sig], gwas_Fin$estim[ind_sig]); abline(0, 1, col = "red")
+plot(gwas_NW$estim[ind_sig], gwas_It$estim[ind_sig]); abline(0, 1, col = "red")
+
+reg <- deming::deming(gwas_NW$estim[ind_sig] ~ gwas_Fin$estim[ind_sig] + 0,
+                      xstd = gwas_Fin$std.err[ind_sig], ystd = gwas_NW$std.err[ind_sig])
+reg
+
+deming::deming(gwas_NW$estim[ind_sig] ~ gwas_It$estim[ind_sig] + 0,
+               xstd = gwas_It$std.err[ind_sig], ystd = gwas_NW$std.err[ind_sig])
+
+
+ind_chr6 <- which(CHR == 6)
+POS2 <- snp_asGeneticPos(CHR[ind_chr6], POS[ind_chr6])
+corr_chr6_NW <- snp_cor(G, ind.row = ind_NW, ind.col = ind_chr6,
+                        infos.pos = POS2, size = 3 / 1000, ncores = NCORES)
+corr_chr6_Fin <- snp_cor(G, ind.row = ind_Fin, ind.col = ind_chr6,
+                         infos.pos = POS2, size = 3 / 1000, ncores = NCORES)
+corr_chr6_It <- snp_cor(G, ind.row = ind_It, ind.col = ind_chr6,
+                        infos.pos = POS2, size = 3 / 1000, ncores = NCORES)
+
+ind <- Matrix::which(corr_chr6_NW^2 > 0.1 & corr_chr6_Fin^2 > 0.1, arr.ind = TRUE)
+ind2 <- ind[ind[, 1] > ind[, 2], ]
+plot(corr_chr6_Fin[ind2], corr_chr6_NW[ind2], pch = 20); abline(0, 1, col = "red")
+
+
+ind3 <- Matrix::which(corr_chr6_NW^2 > 0.1 & corr_chr6_It^2 > 0.1, arr.ind = TRUE)
+ind4 <- ind3[ind3[, 1] > ind3[, 2], ]
+plot(corr_chr6_It[ind4], corr_chr6_NW[ind4], pch = 20); abline(0, 1, col = "red")
+
+plot(corr_chr6_It[ind4], corr_chr6_Fin[ind4], pch = 20); abline(0, 1, col = "red")
+
+
+ind_UK2 <- which(ldist_to_centers[, "UK2"] < 5 & pop2 == "UK2")
+
+gwas_UK1 <- run_GWAS(ind_UK1)
+gwas_UK2 <- run_GWAS(ind_UK2)
+
+plot_grid(
+  snp_manhattan(gwas_UK1,  CHR, POS, npoints = 20e3) +
+    ggplot2::scale_y_log10(),
+  snp_manhattan(gwas_UK2, CHR, POS, npoints = 20e3) +
+    ggplot2::scale_y_log10(),
+  ncol = 1
+)
+corr_chr6_UK2 <- snp_cor(G, ind.row = ind_UK2, ind.col = ind_chr6,
+                         infos.pos = POS2, size = 3 / 1000, ncores = NCORES)
+plot(corr_chr6_UK1[ind2], corr_chr6_UK2[ind2], pch = 20); abline(0, 1, col = "red")
+
+
+ind_Net <- which(ldist_to_centers[, "Netherlands"] < 5 & pop2 == "Netherlands")
+corr_chr6_Net <- snp_cor(G, ind.row = ind_Net, ind.col = ind_chr6,
+                         infos.pos = POS2, size = 3 / 1000, ncores = NCORES)
+plot(corr_chr6_Net[ind2], corr_chr6_UK2[ind2], pch = 20); abline(0, 1, col = "red")
+
+
+ind_UK1 <- sample(which(ldist_to_centers[, "UK1"] < 5 & pop2 == "UK1"), length(ind_Net))
+corr_chr6_UK1 <- snp_cor(G, ind.row = ind_UK1, ind.col = ind_chr6,
+                         infos.pos = POS2, size = 3 / 1000, ncores = NCORES)
+plot(corr_chr6_UK1[ind2], corr_chr6_UK2[ind2], pch = 20); abline(0, 1, col = "red")
+
+
+ind_max <- which.min(predict(gwas_NW))
+ind5 <- ind_max + c(-4:2, 4, 10)
+plot_grid(
+  snp_manhattan(gwas_NW[ind5, ],  CHR[ind5], POS[ind5]) +
+    ggplot2::scale_y_log10(),
+  snp_manhattan(gwas_Fin[ind5, ], CHR[ind5], POS[ind5]) +
+    ggplot2::scale_y_log10(),
+  snp_manhattan(gwas_It[ind5, ],  CHR[ind5], POS[ind5]) +
+    ggplot2::scale_y_log10(),
+  ncol = 1
+)
+
+ind6 <- match(ind5, ind_chr6)
+options(max.print = 200)
+rbind(
+  corr_chr6_NW[ind6, ind6[5]],
+  corr_chr6_Fin[ind6, ind6[5]],
+  corr_chr6_It[ind6, ind6[5]])
diff --git a/tmp-tests/test-linear-solver2.R b/tmp-tests/test-linear-solver2.R
new file mode 100644
index 0000000..fd23dbd
--- /dev/null
+++ b/tmp-tests/test-linear-solver2.R
@@ -0,0 +1,196 @@
+library(bigsnpr)
+
+chr22 <- snp_attach("../Dubois2010_data/celiac_chr22.rds")
+G <- chr22$genotypes$copy(code = c(0, 1, 2, 0, rep(NA, 252)))
+dim(G) #  11402 x 4945
+big_counts(G, ind.col = 1:10)
+CHR <- chr22$map$chromosome
+POS <- chr22$map$physical.pos
+stats <- big_scale()(G)
+
+POS2 <- snp_asGeneticPos(CHR, POS, dir = "tmp-data/")
+plot(POS, POS2, pch = 20)
+
+corr <- runonce::save_run(
+  snp_cor(chr22$genotypes, infos.pos = POS2, size = 3 / 1000, ncores = 6),
+  file = "tmp-data/corr_chr22.rds"
+)
+
+# hist(log10(abs(corr@x)))
+
+
+
+# Simu phenotype
+set.seed(1)
+(M <- round(ncol(G) * 10^-runif(1, 1, 2)))  # 268
+simu <- snp_simuPheno(G, h2 = 0.2, M = M)
+
+y <- simu$pheno
+
+# GWAS
+set.seed(1)
+ind.gwas <- sample(nrow(G), 8e3)
+ind.val <- setdiff(rows_along(G), ind.gwas)
+gwas <- big_univLinReg(G, y[ind.gwas], ind.train = ind.gwas)
+plot(gwas)
+plot(gwas, type = "Manhattan")
+
+df_beta <- data.frame(beta = gwas$estim, beta_se = gwas$std.err,
+                      n_eff = length(ind.gwas))
+
+
+
+# input parameters
+n_vec <- df_beta$n_eff
+beta_hat <- with(df_beta, beta / sqrt(n_eff * beta_se^2 + beta^2))
+
+mean_ld <- mean(Matrix::colSums(corr^2))
+
+burn_in <- 100
+num_iter <- 100
+
+p_init <- 0.01
+h2_init <- 0.2
+
+
+#### LDpred2-auto ####
+
+corr2 <- as.matrix(corr)
+# corr2 <- cov2cor(glasso$w)
+
+m <- length(beta_hat)
+
+{
+  curr_beta <- dotprods <- avg_beta <- avg_postp <- avg_beta_hat <-
+    rep(0, m)
+  num_iter_tot <- burn_in + num_iter
+  p_est <- h2_est <- rep(NA_real_, num_iter_tot)
+
+  id <- matrix(0, m, 1)
+
+  cur_h2_est <- 0
+  p <- p_init
+  h2 <- h2_init
+  gap0 <- crossprod(beta_hat)
+  shrink_corr <- 0.99
+
+  for (k in seq_len(num_iter_tot)) {
+
+    inv_odd_p = (1 - p) / p
+    sigma2 = h2 / (m * p)
+    gap = 0
+    nb_causal <- 0
+
+    for (j in seq_len(m)) {
+
+      # print(j)
+
+      dotprod = dotprods[j];
+      resid = beta_hat[j] - dotprod;
+      gap =  gap + resid * resid;
+      res_beta_hat_j = beta_hat[j] + shrink_corr * (curr_beta[j] - dotprod);
+
+      C1 = sigma2 * n_vec[j];
+      C2 = 1 / (1 + 1 / C1);
+      C3 = C2 * res_beta_hat_j;
+      C4 = C2 / n_vec[j];
+
+      postp = 1 / (1 + inv_odd_p * sqrt(1 + C1) * exp(-C3 * C3 / C4 / 2));
+
+      dotprod_shrunk = shrink_corr * dotprod + (1 - shrink_corr) * curr_beta[j];
+
+      if (k > burn_in) {
+        avg_postp[j]    = avg_postp[j] + postp;
+        avg_beta[j]     = avg_beta[j] + C3 * postp;
+        avg_beta_hat[j] = avg_beta_hat[j] + dotprod_shrunk;
+      }
+
+      if (postp > runif(1)) {
+
+        samp_beta = rnorm(1, C3, sqrt(C4));
+
+        diff = samp_beta - curr_beta[j];
+        curr_beta[j] = samp_beta;
+        nb_causal <- nb_causal + 1
+
+      } else {
+        diff = -curr_beta[j];
+        curr_beta[j] = 0;
+      }
+
+      if (diff != 0) {
+        cur_h2_est = cur_h2_est + diff * (2 * dotprod_shrunk + diff);
+        dotprods = dotprods + corr2[, j] * diff
+        # id[[j]] <- diff
+        # dotprods = dotprods + Matrix::solve(prec2, id)[, 1]
+        # id[[j]] <- 0
+      }
+    }
+
+    if (gap > gap0) stop("Divergence!")
+
+    p = rbeta(1, 1 + nb_causal / mean_ld, 1 + (m - nb_causal) / mean_ld)
+    h2 = cur_h2_est
+    # print(shrink_corr <- crossprod(dotprods, beta_hat) / crossprod(dotprods))
+    # if (shrink_corr > 0.99) shrink_corr <- 0.99
+    # if (shrink_corr < 0.2)  shrink_corr <- 0.2
+
+
+    print(c(k, p, h2))
+    p_est[k]  = p;
+    h2_est[k] = h2;
+  }
+}
+
+
+plot(p_est, log = "y", pch = 20); abline(h = M / m, col = "red")
+plot(h2_est, pch = 20); abline(h = 0.2, col = "red")
+
+pred <- big_prodVec(G, avg_beta, ind.row = ind.val)
+cor(pred, y[ind.val])^2 # 0.143 / 0.143
+# using the normal corr leads to divergence
+
+
+
+id <- rep(0, m)
+id[[100]] <- 1
+id2 <- as.matrix(id)
+
+corr3 <- Matrix::solve(prec)
+prec2 <- as(prec, "dgCMatrix")
+prec3 <- as_SFBM(prec)
+chol_prec <- Matrix::Cholesky(prec2)
+object.size(prec) / object.size(chol_prec) # 2.4
+
+chol_prec2 <- SparseM::chol(prec2)
+object.size(prec) / object.size(chol_prec2) # 1.7
+
+microbenchmark::microbenchmark(
+  Matrix::solve(prec, id)[, 1],
+  Matrix::solve(prec, id2)[, 1],
+  # Matrix::solve(prec2, id)[, 1],
+  Matrix::solve(prec2, id2)[, 1],
+  Matrix::solve(chol_prec, id2)[, 1],
+  bigsparser::sp_solve_sym(prec3, id),
+  times = 10,
+  check = "equal"
+)
+# Unit: milliseconds
+# expr                                     min      lq     mean   median      uq     max neval
+# Matrix::solve(prec, id)[, 1]         49.4297 49.7374 52.13781 51.57370 54.8071 55.1914    10
+# Matrix::solve(prec, id2)[, 1]        48.8244 49.8240 50.76866 50.41285 51.0342 54.5098    10
+# Matrix::solve(prec2, id2)[, 1]       14.4714 14.6538 15.56718 14.73695 16.3233 18.6680    10
+# Matrix::solve(chol_prec, id2)[, 1]   25.9950 26.9176 27.67761 27.22135 28.8990 29.7607    10
+# bigsparser::sp_solve_sym(prec3, id)  51.0848 52.2200 57.17864 53.94215 61.4573 72.8552    10
+
+
+microbenchmark::microbenchmark(
+  Matrix::solve(prec2, id2)[, 1],
+  SparseM::backsolve(chol_prec2, SparseM::forwardsolve(chol_prec2, id), twice = FALSE),
+  SparseM::backsolve(chol_prec2, id),
+  # check = "equal",
+  times = 10
+)
+
+prec3 <- as.matrix(prec2)
+solve_least_squares_qr(prec3, id2)
diff --git a/tmp-tests/test-postp-multipop.R b/tmp-tests/test-postp-multipop.R
new file mode 100644
index 0000000..a0dde42
--- /dev/null
+++ b/tmp-tests/test-postp-multipop.R
@@ -0,0 +1,141 @@
+library(bigsnpr)
+library(dplyr)
+
+celiac <- snp_attach("../Dubois2010_data/FinnuncorrNLITUK1UK3hap300_QC_norel.rds")
+dim(G <- celiac$genotypes)  # 15155 x 281122
+CHR <- celiac$map$chromosome
+POS <- celiac$map$physical.pos
+y <- celiac$fam$affection - 1
+sex <- celiac$fam$sex
+NCORES <- nb_cores()
+ind_chr6 <- which(CHR == 6)
+
+pop <- snp_getSampleInfos(celiac, "../Dubois2010_data/FinnNLITUK1UK3.clusterv2")[[1]]
+pop2 <- c("Netherlands", "Italy", "UK1", "UK2", "Finland")[pop]
+table(pop2, exclude = NULL)
+# Finland       Italy Netherlands         UK1         UK2
+#    2436        1035        1623        3325        6736
+
+big_counts(G, ind.col = 1:10)
+G2 <- snp_fastImputeSimple(G, method = "mean2", ncores = NCORES)
+
+obj.svd <- runonce::save_run(
+  snp_autoSVD(G2, CHR, POS, ncores = NCORES),
+  file = "tmp-data/celiac_svd.rds"
+)
+plot(obj.svd)
+plot(obj.svd, type = "scores", scores = 1:6, coef = 0.6)
+PC <- predict(obj.svd)[, 1:6]
+
+pop_centers <- bigutilsr::geometric_median(PC, by_grp = pop2)
+
+ldist_to_centers <- apply(pop_centers, 1, function(center) {
+  log(rowSums(sweep(PC, 2, center, '-')^2))
+})
+
+ind_Fin <- which(ldist_to_centers[, "Finland"] < 7 & pop2 == "Finland")
+ind_It <- which(ldist_to_centers[, "Italy"] < 6 & pop2 == "Italy")
+
+corr_It <- runonce::save_run(file = "tmp-data/LD_It_chr6.rds", {
+  POS2 <- snp_asGeneticPos(CHR[ind_chr6], POS[ind_chr6])
+  corr0 <- snp_cor(G, ind.row = ind_It, ind.col = ind_chr6,
+                   size = 3 / 1000, infos.pos = POS2, ncores = NCORES)
+  as_SFBM(corr0, backingfile = "tmp-data/LD_It_chr6", compact = TRUE)
+})
+
+corr_Fin <- runonce::save_run(file = "tmp-data/LD_Fin_chr6.rds", {
+  POS2 <- snp_asGeneticPos(CHR[ind_chr6], POS[ind_chr6])
+  corr0 <- snp_cor(G, ind.row = ind_Fin, ind.col = ind_chr6,
+                   size = 3 / 1000, infos.pos = POS2, ncores = NCORES)
+  as_SFBM(corr0, backingfile = "tmp-data/LD_Fin_chr6", compact = TRUE)
+})
+
+ind_ItFin <- c(ind_It, ind_Fin)
+corr_ItFin <- runonce::save_run(file = "tmp-data/LD_ItFin_chr6.rds", {
+  POS2 <- snp_asGeneticPos(CHR[ind_chr6], POS[ind_chr6])
+  corr0 <- snp_cor(G, ind.row = ind_ItFin, ind.col = ind_chr6,
+                   size = 3 / 1000, infos.pos = POS2, ncores = NCORES)
+  as_SFBM(corr0, backingfile = "tmp-data/LD_ItFin_chr6", compact = TRUE)
+})
+
+get_postp <- function(pop) {
+  ind <- get(paste0("ind_", pop))
+  gwas <- big_univLogReg(G2, y[ind], ind.train = ind, ind.col = ind_chr6,
+                         covar.train = cbind(sex, PC)[ind, ])
+  df_beta <- dplyr::transmute(gwas, beta = estim, beta_se = std.err, n_eff = length(ind))
+
+  corr <- get(paste0("corr_", pop))
+  ldpred_auto0 <- snp_ldpred2_auto(corr, df_beta, h2_init = 0.2, vec_p_init = 0.01,
+                                   burn_in = 100, num_iter = 200, verbose = TRUE,
+                                   allow_jump_sign = FALSE, shrink_corr = 0.95,
+                                   use_MLE = FALSE)[[1]]
+  ldpred_auto0$postp_est
+}
+
+postp_It <- get_postp("It")
+summary(postp_It)
+#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
+# 0.003458 0.003690 0.004188 0.006615 0.005710 0.897038
+
+postp_Fin <- get_postp("Fin")
+summary(postp_Fin)
+#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
+# 0.006014 0.006613 0.007705 0.013554 0.010911 1.000000
+
+postp_ItFin <- get_postp("ItFin")
+summary(postp_ItFin)
+#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
+# 0.005426 0.005999 0.007040 0.014546 0.010238 1.000000
+
+
+library(ggplot2)
+library(dplyr)
+tibble(postp_It, postp_Fin, postp_ItFin) %>%
+  tidyr::pivot_longer(cols = 1:3) %>%
+  ggplot(aes(value, fill = name)) +
+  geom_density(alpha = 0.4) +
+  scale_x_log10()
+
+
+get_beta_hat <- function(ind) {
+  gwas <- big_univLogReg(G2, y[ind], ind.train = ind, ind.col = ind_chr6,
+                         covar.train = cbind(sex, PC)[ind, ])
+  with(gwas, estim / sqrt(length(ind) * std.err^2 + estim^2))
+}
+beta_hat <- cbind(get_beta_hat(ind_It), get_beta_hat(ind_Fin))
+plot(beta_hat)
+
+var <- cbind(big_colstats(G2, ind.row = ind_It, ind.col = ind_chr6)$var,
+             big_colstats(G2, ind.row = ind_Fin, ind.col = ind_chr6)$var)
+plot(var)
+
+mean_ld <- mean(
+  bigsnpr:::ld_scores_sfbm(corr_It, ind_sub = cols_along(corr_It) - 1L, ncores = 2))
+
+Rcpp::sourceCpp("tmp-tests/proto-ldpred2x.cpp")
+
+ldpred_auto <- ldpred2x_gibbs(
+  list_corr    = list(corr_It, corr_Fin),
+  beta_hat     = beta_hat,
+  n_vec        = matrix(c(length(ind_It), length(ind_Fin)),
+                        nrow = ncol(corr_It), ncol = 2, byrow = TRUE),
+  log_var      = log(var),
+  ind_sub      = cols_along(corr_It) - 1L,
+  p_init       = 0.01,
+  h2_init      = 0.5,
+  burn_in      = 100,
+  num_iter     = 200,
+  verbose      = TRUE,
+  no_jump_sign = TRUE,
+  shrink_corr  = 0.95,
+  use_mle      = FALSE,
+  alpha_bounds = c(-1.5, 0.5) + 1,
+  mean_ld      = mean_ld
+)
+
+postp_x <- ldpred_auto$postp_est
+
+plot(postp_x, postp_ItFin, pch = 20); abline(0, 1, col = "red", lwd = 2)
+plot(tibble(postp_It, postp_Fin, postp_ItFin, postp_x))
+plot(postp_x, postp_It, pch = 20); abline(0, 1, col = "red", lwd = 2)
+plot(postp_x, postp_Fin, pch = 20); abline(0, 1, col = "red", lwd = 2)