From 2bcc457bf7ee7f8e4256579688d57663b67a3b74 Mon Sep 17 00:00:00 2001 From: Freguglia Date: Thu, 28 May 2026 19:34:17 -0300 Subject: [PATCH 1/4] Refactor conditional probability functions to use log-space calculations and improve numerical stability --- DESCRIPTION | 2 +- src/conditional_probability.cpp | 14 ++++-- src/pseudo_likelihood.cpp | 81 +++++++++++++++++++++++++++------ 3 files changed, 78 insertions(+), 19 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 35aa12c..b0391c5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -29,7 +29,7 @@ LinkingTo: Rcpp, RcppArmadillo Roxygen: list(markdown = TRUE) RdMacros: Rdpack -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.3 Suggests: testthat (>= 2.1.0), covr, diff --git a/src/conditional_probability.cpp b/src/conditional_probability.cpp index e6d0efd..80a0f37 100644 --- a/src/conditional_probability.cpp +++ b/src/conditional_probability.cpp @@ -15,7 +15,7 @@ NumericVector conditional_probabilities_mrf(const IntegerMatrix &Z, const int N, const int M, const int n_R, const int C){ - NumericVector probs(C+1); + NumericVector log_pots(C+1); double this_prob; int dx, dy; int x = position[0] -1; int y = position[1] -1; @@ -29,8 +29,11 @@ NumericVector conditional_probabilities_mrf(const IntegerMatrix &Z, if(0 <= x-dx && x-dx < N && 0 <= y-dy && y-dy < M){ this_prob = this_prob + theta(Z(x-dx, y-dy), value, i);} } - probs[value] = exp(this_prob); + log_pots[value] = this_prob; } + // log-sum-exp trick: subtract max before exponentiation to prevent overflow + double max_lp = max(log_pots); + NumericVector probs = exp(log_pots - max_lp); return(probs/sum(probs)); } @@ -44,7 +47,7 @@ NumericVector conditional_probabilities_mrf_sub(const IntegerMatrix &Z, const int N, const int M, const int n_R, const int C){ - NumericVector probs(C+1); + NumericVector log_pots(C+1); double this_prob; int dx, dy; int x = position[0] -1; int y = position[1] -1; @@ -60,8 +63,11 @@ NumericVector conditional_probabilities_mrf_sub(const IntegerMatrix &Z, if(sub_mat(x-dx, y-dy)){ this_prob = this_prob + theta(Z(x-dx, y-dy), value, i);}} } - probs[value] = exp(this_prob); + log_pots[value] = this_prob; } + // log-sum-exp trick: subtract max before exponentiation to prevent overflow + double max_lp = max(log_pots); + NumericVector probs = exp(log_pots - max_lp); return(probs/sum(probs)); } diff --git a/src/pseudo_likelihood.cpp b/src/pseudo_likelihood.cpp index 105c8dc..4129d88 100644 --- a/src/pseudo_likelihood.cpp +++ b/src/pseudo_likelihood.cpp @@ -3,6 +3,69 @@ using namespace Rcpp; +// Computes log P(Z_ij = zij | Z_rest) directly in log-space: +// log P = s_zij - log_sum_exp(s_0, ..., s_C) +// where s_v = sum of theta potentials for value v with all neighbors. +// Avoids: per-pixel NumericVector allocation, redundant exp/log roundtrip, +// and numerical overflow via the log-sum-exp trick. +static inline double log_cond_pixel(const IntegerMatrix &Z, + const IntegerMatrix &R, + const arma::fcube &theta, + const int N, const int M, + const int n_R, const int C, + const int i, const int j, + std::vector &log_probs){ + int dx, dy; + for(int v = 0; v <= C; v++){ + double lp = 0.0; + for(int r = 0; r < n_R; r++){ + dx = R(r, 0); dy = R(r, 1); + if(i+dx >= 0 && i+dx < N && j+dy >= 0 && j+dy < M) + lp += theta(v, Z(i+dx, j+dy), r); + if(i-dx >= 0 && i-dx < N && j-dy >= 0 && j-dy < M) + lp += theta(Z(i-dx, j-dy), v, r); + } + log_probs[v] = lp; + } + + // log-sum-exp trick: subtract max before exponentiation to prevent overflow + double max_lp = *std::max_element(log_probs.begin(), log_probs.end()); + double sum_exp = 0.0; + for(int v = 0; v <= C; v++) + sum_exp += std::exp(log_probs[v] - max_lp); + + return log_probs[Z(i, j)] - max_lp - std::log(sum_exp); +} + +static inline double log_cond_pixel_sub(const IntegerMatrix &Z, + const LogicalMatrix &sub_mat, + const IntegerMatrix &R, + const arma::fcube &theta, + const int N, const int M, + const int n_R, const int C, + const int i, const int j, + std::vector &log_probs){ + int dx, dy; + for(int v = 0; v <= C; v++){ + double lp = 0.0; + for(int r = 0; r < n_R; r++){ + dx = R(r, 0); dy = R(r, 1); + if(i+dx >= 0 && i+dx < N && j+dy >= 0 && j+dy < M && sub_mat(i+dx, j+dy)) + lp += theta(v, Z(i+dx, j+dy), r); + if(i-dx >= 0 && i-dx < N && j-dy >= 0 && j-dy < M && sub_mat(i-dx, j-dy)) + lp += theta(Z(i-dx, j-dy), v, r); + } + log_probs[v] = lp; + } + + double max_lp = *std::max_element(log_probs.begin(), log_probs.end()); + double sum_exp = 0.0; + for(int v = 0; v <= C; v++) + sum_exp += std::exp(log_probs[v] - max_lp); + + return log_probs[Z(i, j)] - max_lp - std::log(sum_exp); +} + // [[Rcpp::depends(RcppArmadillo)]] // [[Rcpp::export]] double log_pl_mrf(const IntegerMatrix Z, @@ -13,16 +76,11 @@ double log_pl_mrf(const IntegerMatrix Z, int n_R = R.nrow(); int C = theta.n_rows - 1; double log_pl = 0.0; - double this_cond_prob = 0.0; - IntegerVector position(2); - int zij = -1; + std::vector log_probs(C + 1); for(int i = 0; i < N; i++){ for(int j = 0; j < M; j++){ - zij = Z(i,j); - position[0] = i+1; position[1] = j+1; - this_cond_prob = conditional_probabilities_mrf(Z, position, R, theta, N, M, n_R, C)(zij); - log_pl += log(this_cond_prob); + log_pl += log_cond_pixel(Z, R, theta, N, M, n_R, C, i, j, log_probs); } } return(log_pl); @@ -39,17 +97,12 @@ double log_pl_mrf_sub(const IntegerMatrix Z, int n_R = R.nrow(); int C = theta.n_rows - 1; double log_pl = 0.0; - double this_cond_prob = 0.0; - IntegerVector position(2); - int zij = -1; + std::vector log_probs(C + 1); for(int i = 0; i < N; i++){ for(int j = 0; j < M; j++){ if(sub_mat(i,j)){ - zij = Z(i,j); - position[0] = i+1; position[1] = j+1; - this_cond_prob = conditional_probabilities_mrf_sub(Z, sub_mat, position, R, theta, N, M, n_R, C)(zij); - log_pl += log(this_cond_prob); + log_pl += log_cond_pixel_sub(Z, sub_mat, R, theta, N, M, n_R, C, i, j, log_probs); } } } From 4b4418c4d67266a137a4e89a34950bd06d16e53c Mon Sep 17 00:00:00 2001 From: Freguglia Date: Thu, 28 May 2026 20:27:13 -0300 Subject: [PATCH 2/4] Bump version to 1.0.2, add OpenMP support in pseudolikelihood calculations, and optimize conditional probability functions for performance --- DESCRIPTION | 2 +- src/Makevars | 2 + src/Makevars.win | 2 + src/conditional_probability.cpp | 52 +++++++++---------- src/pseudo_likelihood.cpp | 88 +++++++++++++++++++++++++-------- 5 files changed, 99 insertions(+), 47 deletions(-) create mode 100644 src/Makevars create mode 100644 src/Makevars.win diff --git a/DESCRIPTION b/DESCRIPTION index b0391c5..48a1301 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: mrf2d Type: Package Title: Markov Random Field Models for Image Analysis -Version: 1.0.1 +Version: 1.0.2 Authors@R: person("Victor", "Freguglia", email = "victorfreguglia@gmail.com", role = c("aut", "cre"), diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..a576f88 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +PKG_CXXFLAGS += $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS += $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..a576f88 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,2 @@ +PKG_CXXFLAGS += $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS += $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/conditional_probability.cpp b/src/conditional_probability.cpp index 80a0f37..54c0549 100644 --- a/src/conditional_probability.cpp +++ b/src/conditional_probability.cpp @@ -15,21 +15,22 @@ NumericVector conditional_probabilities_mrf(const IntegerMatrix &Z, const int N, const int M, const int n_R, const int C){ - NumericVector log_pots(C+1); - double this_prob; + NumericVector log_pots(C+1, 0.0); int dx, dy; int x = position[0] -1; int y = position[1] -1; - for(int value = 0; value <= C; value++){ - this_prob = 0.0; - for(int i = 0; i < n_R; i++){ - dx = R(i,0); dy = R(i,1); - if(0 <= x+dx && x+dx < N && 0 <= y+dy && y+dy < M){ - this_prob = this_prob + theta(value, Z(x+dx, y+dy), i);} - if(0 <= x-dx && x-dx < N && 0 <= y-dy && y-dy < M){ - this_prob = this_prob + theta(Z(x-dx, y-dy), value, i);} + for(int i = 0; i < n_R; i++){ + dx = R(i,0); dy = R(i,1); + if(0 <= x+dx && x+dx < N && 0 <= y+dy && y+dy < M){ + int nb = Z(x+dx, y+dy); + for(int value = 0; value <= C; value++) + log_pots[value] += theta(value, nb, i); + } + if(0 <= x-dx && x-dx < N && 0 <= y-dy && y-dy < M){ + int nb = Z(x-dx, y-dy); + for(int value = 0; value <= C; value++) + log_pots[value] += theta(nb, value, i); } - log_pots[value] = this_prob; } // log-sum-exp trick: subtract max before exponentiation to prevent overflow double max_lp = max(log_pots); @@ -47,23 +48,24 @@ NumericVector conditional_probabilities_mrf_sub(const IntegerMatrix &Z, const int N, const int M, const int n_R, const int C){ - NumericVector log_pots(C+1); - double this_prob; + NumericVector log_pots(C+1, 0.0); int dx, dy; int x = position[0] -1; int y = position[1] -1; - for(int value = 0; value <= C; value++){ - this_prob = 0.0; - for(int i = 0; i < n_R; i++){ - dx = R(i,0); dy = R(i,1); - if(0 <= x+dx && x+dx < N && 0 <= y+dy && y+dy < M){ - if(sub_mat(x+dx, y+dy)){ - this_prob = this_prob + theta(value, Z(x+dx, y+dy), i);}} - if(0 <= x-dx && x-dx < N && 0 <= y-dy && y-dy < M){ - if(sub_mat(x-dx, y-dy)){ - this_prob = this_prob + theta(Z(x-dx, y-dy), value, i);}} - } - log_pots[value] = this_prob; + for(int i = 0; i < n_R; i++){ + dx = R(i,0); dy = R(i,1); + if(0 <= x+dx && x+dx < N && 0 <= y+dy && y+dy < M){ + if(sub_mat(x+dx, y+dy)){ + int nb = Z(x+dx, y+dy); + for(int value = 0; value <= C; value++) + log_pots[value] += theta(value, nb, i); + }} + if(0 <= x-dx && x-dx < N && 0 <= y-dy && y-dy < M){ + if(sub_mat(x-dx, y-dy)){ + int nb = Z(x-dx, y-dy); + for(int value = 0; value <= C; value++) + log_pots[value] += theta(nb, value, i); + }} } // log-sum-exp trick: subtract max before exponentiation to prevent overflow double max_lp = max(log_pots); diff --git a/src/pseudo_likelihood.cpp b/src/pseudo_likelihood.cpp index 4129d88..1b402cc 100644 --- a/src/pseudo_likelihood.cpp +++ b/src/pseudo_likelihood.cpp @@ -1,5 +1,8 @@ #include #include "mrf2d.h" +#ifdef _OPENMP +#include +#endif using namespace Rcpp; @@ -16,19 +19,27 @@ static inline double log_cond_pixel(const IntegerMatrix &Z, const int i, const int j, std::vector &log_probs){ int dx, dy; - for(int v = 0; v <= C; v++){ - double lp = 0.0; - for(int r = 0; r < n_R; r++){ - dx = R(r, 0); dy = R(r, 1); - if(i+dx >= 0 && i+dx < N && j+dy >= 0 && j+dy < M) - lp += theta(v, Z(i+dx, j+dy), r); - if(i-dx >= 0 && i-dx < N && j-dy >= 0 && j-dy < M) - lp += theta(Z(i-dx, j-dy), v, r); + + for(int v = 0; v <= C; v++) log_probs[v] = 0.0; + + // Neighbors outer, values inner: + // - Z(neighbor) read once per direction instead of C+1 times + // - bounds check done once per direction + // - theta(v, nb, r) with v sequential = cache-friendly (Armadillo col-major) + for(int r = 0; r < n_R; r++){ + dx = R(r, 0); dy = R(r, 1); + if(i+dx >= 0 && i+dx < N && j+dy >= 0 && j+dy < M){ + int nb = Z(i+dx, j+dy); + for(int v = 0; v <= C; v++) + log_probs[v] += theta(v, nb, r); + } + if(i-dx >= 0 && i-dx < N && j-dy >= 0 && j-dy < M){ + int nb = Z(i-dx, j-dy); + for(int v = 0; v <= C; v++) + log_probs[v] += theta(nb, v, r); } - log_probs[v] = lp; } - // log-sum-exp trick: subtract max before exponentiation to prevent overflow double max_lp = *std::max_element(log_probs.begin(), log_probs.end()); double sum_exp = 0.0; for(int v = 0; v <= C; v++) @@ -46,16 +57,21 @@ static inline double log_cond_pixel_sub(const IntegerMatrix &Z, const int i, const int j, std::vector &log_probs){ int dx, dy; - for(int v = 0; v <= C; v++){ - double lp = 0.0; - for(int r = 0; r < n_R; r++){ - dx = R(r, 0); dy = R(r, 1); - if(i+dx >= 0 && i+dx < N && j+dy >= 0 && j+dy < M && sub_mat(i+dx, j+dy)) - lp += theta(v, Z(i+dx, j+dy), r); - if(i-dx >= 0 && i-dx < N && j-dy >= 0 && j-dy < M && sub_mat(i-dx, j-dy)) - lp += theta(Z(i-dx, j-dy), v, r); + + for(int v = 0; v <= C; v++) log_probs[v] = 0.0; + + for(int r = 0; r < n_R; r++){ + dx = R(r, 0); dy = R(r, 1); + if(i+dx >= 0 && i+dx < N && j+dy >= 0 && j+dy < M && sub_mat(i+dx, j+dy)){ + int nb = Z(i+dx, j+dy); + for(int v = 0; v <= C; v++) + log_probs[v] += theta(v, nb, r); + } + if(i-dx >= 0 && i-dx < N && j-dy >= 0 && j-dy < M && sub_mat(i-dx, j-dy)){ + int nb = Z(i-dx, j-dy); + for(int v = 0; v <= C; v++) + log_probs[v] += theta(nb, v, r); } - log_probs[v] = lp; } double max_lp = *std::max_element(log_probs.begin(), log_probs.end()); @@ -67,6 +83,7 @@ static inline double log_cond_pixel_sub(const IntegerMatrix &Z, } // [[Rcpp::depends(RcppArmadillo)]] +// [[Rcpp::plugins(openmp)]] // [[Rcpp::export]] double log_pl_mrf(const IntegerMatrix Z, const IntegerMatrix R, @@ -76,17 +93,31 @@ double log_pl_mrf(const IntegerMatrix Z, int n_R = R.nrow(); int C = theta.n_rows - 1; double log_pl = 0.0; - std::vector log_probs(C + 1); +#ifdef _OPENMP + #pragma omp parallel reduction(+:log_pl) + { + std::vector log_probs(C + 1); // thread-local buffer + #pragma omp for collapse(2) schedule(static) + for(int i = 0; i < N; i++){ + for(int j = 0; j < M; j++){ + log_pl += log_cond_pixel(Z, R, theta, N, M, n_R, C, i, j, log_probs); + } + } + } +#else + std::vector log_probs(C + 1); for(int i = 0; i < N; i++){ for(int j = 0; j < M; j++){ log_pl += log_cond_pixel(Z, R, theta, N, M, n_R, C, i, j, log_probs); } } +#endif return(log_pl); } // [[Rcpp::depends(RcppArmadillo)]] +// [[Rcpp::plugins(openmp)]] // [[Rcpp::export]] double log_pl_mrf_sub(const IntegerMatrix Z, const LogicalMatrix sub_mat, @@ -97,8 +128,22 @@ double log_pl_mrf_sub(const IntegerMatrix Z, int n_R = R.nrow(); int C = theta.n_rows - 1; double log_pl = 0.0; - std::vector log_probs(C + 1); +#ifdef _OPENMP + #pragma omp parallel reduction(+:log_pl) + { + std::vector log_probs(C + 1); // thread-local buffer + #pragma omp for collapse(2) schedule(static) + for(int i = 0; i < N; i++){ + for(int j = 0; j < M; j++){ + if(sub_mat(i,j)){ + log_pl += log_cond_pixel_sub(Z, sub_mat, R, theta, N, M, n_R, C, i, j, log_probs); + } + } + } + } +#else + std::vector log_probs(C + 1); for(int i = 0; i < N; i++){ for(int j = 0; j < M; j++){ if(sub_mat(i,j)){ @@ -106,5 +151,6 @@ double log_pl_mrf_sub(const IntegerMatrix Z, } } } +#endif return(log_pl); } From 3047bfadc912dc7d0c4a41b2f211474321e5e586 Mon Sep 17 00:00:00 2001 From: Freguglia Date: Thu, 28 May 2026 20:33:45 -0300 Subject: [PATCH 3/4] Fix Makevars and Makevars.win to use assignment instead of append for OpenMP flags --- src/Makevars | 4 ++-- src/Makevars.win | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Makevars b/src/Makevars index a576f88..139b1d1 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,2 +1,2 @@ -PKG_CXXFLAGS += $(SHLIB_OPENMP_CXXFLAGS) -PKG_LIBS += $(SHLIB_OPENMP_CXXFLAGS) +PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/Makevars.win b/src/Makevars.win index a576f88..139b1d1 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,2 +1,2 @@ -PKG_CXXFLAGS += $(SHLIB_OPENMP_CXXFLAGS) -PKG_LIBS += $(SHLIB_OPENMP_CXXFLAGS) +PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) From 8de6f66f90930673c0b6443d8d7f6c7198a67a26 Mon Sep 17 00:00:00 2001 From: Freguglia Date: Thu, 28 May 2026 20:40:28 -0300 Subject: [PATCH 4/4] Update upload-artifact action to v4 in test coverage workflow --- .github/workflows/test-coverage.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index d31e1a7..300d6c0 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -41,7 +41,7 @@ jobs: - name: Upload test results if: failure() - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: coverage-test-failures path: ${{ runner.temp }}/package