Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BayesSpace v1.11.3 #111

Open
wants to merge 39 commits into
base: master
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
39 commits
Select commit Hold shift + click to select a range
6b73168
fix: bugs when reading H5
Jun 3, 2023
3c8d3a9
fix: bugs
Jun 4, 2023
5f73f19
fix: bugs
Jun 5, 2023
86de359
fix: added sinity check and core tuning
senbaikang Jul 5, 2023
3c0f502
doc: updated doc
senbaikang Jul 5, 2023
40fdaa1
fix: fixed bugs & faster code
senbaikang Sep 7, 2023
e20e5a3
feat: adaptive mcmc
senbaikang Sep 11, 2023
d6abb37
feat: adaptive mcmc
senbaikang Sep 20, 2023
7186e72
feat: adaptive MCMC
senbaikang Sep 22, 2023
94a9df0
feat: label adjustment
senbaikang Sep 25, 2023
7129c45
fix: bugs
senbaikang Oct 4, 2023
e25d9b2
feat: option for restricted adaptive MCMC
senbaikang Oct 5, 2023
a9c094e
fix: [v1.11.7] bugs in func
senbaikang Dec 8, 2023
b5bd8f2
fix: bugs
senbaikang Feb 6, 2024
f69e023
fix: [v1.11.8] bugs in processing z
senbaikang Feb 16, 2024
6d02fdf
fix: [v1.11.8] bugs in processing z
senbaikang Feb 16, 2024
483645e
fix: bug
senbaikang Feb 19, 2024
65736e4
fix: updated outdated functions
senbaikang Feb 19, 2024
98c1d15
fix: updated outdated functions
senbaikang Feb 19, 2024
8d34497
fix: updated outdated functions
senbaikang Feb 19, 2024
2190e56
feat: [v1.11.10] refactored neighbor finding for subspots
senbaikang Mar 7, 2024
554b98b
fix: bugs
senbaikang Mar 7, 2024
79e11b9
Support for VisiumHD (#1)
senbaikang Mar 7, 2024
fc0d4c7
fix: bugs
senbaikang Mar 8, 2024
4a347ce
fix: cleaned up coords assignment for subspots
senbaikang Mar 13, 2024
6140799
fix: bug in plotting clusters
senbaikang Mar 13, 2024
06b5cf6
feat: [v1.13.1] plot clusters align with the image
senbaikang Mar 14, 2024
ec91302
fix: bugs
senbaikang Mar 14, 2024
d7ad525
feat: [v1.13.2] parallelized PCA
senbaikang Mar 18, 2024
6b28b43
feat: [v1.13.2] removed unrelated probes; parallized modeling of gene…
senbaikang Mar 19, 2024
db2d58d
feat: [v1.13.4] thining in c++ code
senbaikang Mar 20, 2024
dff61c8
fix: bugs
senbaikang Mar 24, 2024
18be765
feat: [v1.13.5] customized number of subspots
senbaikang Mar 25, 2024
9a5f742
fix: bugs
senbaikang Mar 25, 2024
d3e9275
doc: updated doc
senbaikang Mar 26, 2024
1fc79ae
fix: typo
senbaikang Apr 2, 2024
b9640f9
feat: [v1.13.6] added a few helper functions
senbaikang Apr 8, 2024
df216c1
fix: bugs
senbaikang Oct 15, 2024
0d6b7ac
doc: update NEWS
senbaikang Oct 15, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
feat: adaptive mcmc
  • Loading branch information
senbaikang committed Sep 20, 2023
commit d6abb37c161ae226ec97391e718f574d310abfcf
84 changes: 53 additions & 31 deletions src/cluster.cpp
Original file line number Diff line number Diff line change
@@ -43,22 +43,24 @@ sig_handler(int _) {
early_stop = 1;
}

double
update_jitter_scale(
double accpet_rate, double target_accept_rate, size_t min_iter,
size_t curr_iter, double curr_jitter_scale
mat
adaptive_mcmc(
double accpetance_rate, double target_acceptance_rate, size_t curr_iter,
const mat &identity_mtx, const mat &adaptive_mtx, const rowvec &samples
) {
if (min_iter >= curr_iter)
return curr_jitter_scale;

const double step_size = std::min(0.01, 1.0 / sqrt(curr_iter));

if (accpet_rate < target_accept_rate)
return curr_jitter_scale - step_size;
else if (accpet_rate > target_accept_rate)
return curr_jitter_scale + step_size;
else
return curr_jitter_scale;
const double step_size =
std::min(1.0, samples.n_elem * std::pow(curr_iter, -2.0 / 3));

const mat sample_mtx = resize(samples, samples.n_elem, 1);

return chol(
adaptive_mtx *
(identity_mtx +
step_size * (accpetance_rate - target_acceptance_rate) * sample_mtx *
sample_mtx.t() / accu(samples % samples)) *
adaptive_mtx.t(),
"lower"
);
}

/* C++ version of the dtrmv BLAS function */
@@ -716,12 +718,19 @@ iterate_deconv(
const IntegerVector qvec = seq_len(q);
const vec zero_vec = zeros<vec>(d);
const vec one_vec = ones<vec>(d);
const mat error_var =
diagmat(one_vec) / d; // d here does not seem to have an effect
const mat error_var = diagmat(one_vec);
jitter_scale /= d;

// For adaptive MCMC
size_t num_accepts = 0;
size_t num_rejects = 0;
std::vector<size_t> num_accepts(n0, 0);
std::vector<size_t> num_rejects(n0, 0);
std::vector<mat> adaptive_mtx(n);
if (jitter_scale == 0.0) {
#pragma omp parallel for
for (int i = 0; i < n; i++) {
adaptive_mtx[i] = diagmat(one_vec);
}
}

// Progress bar
indicators::show_console_cursor(false);
@@ -747,9 +756,9 @@ iterate_deconv(

#pragma omp parallel shared( \
early_stop, Y, Y_new, error, error_var, num_accepts, num_rejects, \
acceptance_prob, j0_vector, subspots, zero_vec, mu_i_long, Y0, \
lambda_i, n0, c, thread_hits, n, z, z_new, neighbors, gamma, mu_i, \
w, log_likelihoods \
adaptive_mtx, acceptance_prob, j0_vector, subspots, zero_vec, \
mu_i_long, Y0, lambda_i, n0, c, thread_hits, n, z, z_new, \
neighbors, gamma, mu_i, w, log_likelihoods \
)
{
for (int i = 1; i < nrep; i++) {
@@ -782,7 +791,8 @@ iterate_deconv(

// Propose new values for Y.
error = rmvnorm(
n, zero_vec, error_var * jitter_scale
n, zero_vec,
jitter_scale == 0.0 ? error_var : error_var * jitter_scale
); // Generate random numbers before entering multithreading.
}

@@ -797,6 +807,15 @@ iterate_deconv(
const mat Y_j_prev = Y.rows(j0_vector * n0 + j0);
mat error_j = error.rows(j0_vector * n0 + j0);

if (jitter_scale == 0.0) {
for (int r = 0; r < subspots; r++) {
error_j.row(r) = trans(
adaptive_mtx[r * n0 + j0] *
resize(error_j.row(r), error_j.n_cols, 1)
);
}
}

// Make sure that the sum of the error terms is zero.
const rowvec error_mean = sum(error_j, 0) / subspots;
for (int r = 0; r < subspots; r++) {
@@ -845,12 +864,21 @@ iterate_deconv(
if (yesUpdate == 1) {
Y.rows(j0_vector * n0 + j0) = Y_new.rows(j0_vector * n0 + j0);

num_accepts++;
num_accepts[j0]++;
updateCounter++;
} else
num_rejects++;
num_rejects[j0]++;

for (int r = 0; r < subspots; r++) {
// Adaptive MCMC.
if (jitter_scale == 0.0 && i > 10)
adaptive_mtx[r * n0 + j0] = adaptive_mcmc(
static_cast<double>(num_accepts[j0]) /
(num_accepts[j0] + num_rejects[j0]),
0.234, i, error_var, adaptive_mtx[r * n0 + j0],
error.row(r * n0 + j0)
);

// Update w.
if (tdist) {
const double w_beta = as_scalar(
@@ -870,12 +898,6 @@ iterate_deconv(
}
}
Ychange[i] = updateCounter * 1.0 / n0;

jitter_scale = update_jitter_scale(
static_cast<double>(num_accepts) / (num_accepts + num_rejects),
0.234, 10, i, jitter_scale
);
jitterScale[i] = jitter_scale;
}

// Update z