diff --git a/.Rbuildignore b/.Rbuildignore
index cc3b90c3..692b5609 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,15 +1,35 @@
-^renv$
-^renv\.lock$
+# RStudio / IDE
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^\.vscode$
 
-^Readme.Rmd$
-^\.github$
-^_pkgdown\.yml$
+# renv
+^renv$
+^renv\.lock$
+
+# pkgdown / docs
 ^docs$
 ^pkgdown$
-^vignettes/introduction_cache
+^_pkgdown\.yml$
+^Readme\.Rmd$
+^vignettes/introduction_cache$
+
+# GitHub / CI
+^\.github$
+
+# R CMD build artifacts
 ^doc$
 ^Meta$
-^\.vscode$
-^dev/
+
+# Development helpers
+^dev$
+
+# ---- C/C++ build artifacts (REQUIRED) ----
+^src/.*\.o$
+^src/.*\.so$
+^src/.*\.dll$
+
+# ---- Generated build files ----
+^src/Makevars$
+^src/Makevars\.win$
+^src/sources\.mk$
diff --git a/.gitignore b/.gitignore
index 65196573..40fc89cc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,6 +5,9 @@
 src/*.o
 src/*.so
 src/*.dll
+src/**/*.o
+src/**/*.so
+src/**/*.dll
 .DS_Store
 /doc/
 /Meta/
diff --git a/R/RcppExports.R b/R/RcppExports.R
index eb462989..520101b2 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -9,18 +9,6 @@ run_bgm_parallel <- function(observations, num_categories, pairwise_scale, edge_
     .Call(`_bgms_run_bgm_parallel`, observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type)
 }
 
-get_explog_switch <- function() {
-    .Call(`_bgms_get_explog_switch`)
-}
-
-rcpp_ieee754_exp <- function(x) {
-    .Call(`_bgms_rcpp_ieee754_exp`, x)
-}
-
-rcpp_ieee754_log <- function(x) {
-    .Call(`_bgms_rcpp_ieee754_log`, x)
-}
-
 sample_omrf_gibbs <- function(no_states, no_variables, no_categories, interactions, thresholds, iter) {
     .Call(`_bgms_sample_omrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, iter)
 }
diff --git a/R/generate_makevars_sources.R b/R/generate_makevars_sources.R
new file mode 100644
index 00000000..46624760
--- /dev/null
+++ b/R/generate_makevars_sources.R
@@ -0,0 +1,26 @@
+cpp <- list.files(
+  "src",
+  pattern = "\\.cpp$",
+  recursive = TRUE,
+  full.names = TRUE
+)
+
+# strip leading "src/"
+cpp <- sub("^src/", "", cpp)
+
+con <- file("src/sources.mk", open = "w")
+
+writeLines(c(
+  "# ------------------------------------------------------------------",
+  "# THIS FILE IS AUTO-GENERATED - DO NOT EDIT",
+  "# Generated by configure",
+  "# To add C++ code, place .cpp files anywhere under src/",
+  "# ------------------------------------------------------------------",
+  "SOURCES = \\"
+), con)
+
+writeLines(paste0("  ", cpp, " \\"), con)
+writeLines("", con)
+
+close(con)
+
diff --git a/configure b/configure
index 2a547070..ede3910e 100755
--- a/configure
+++ b/configure
@@ -1,9 +1,12 @@
 #!/bin/sh
 
-# Get flags from RcppParallel
+# RcppParallel flags
 RCPP_PARALLEL_CPPFLAGS=`"${R_HOME}/bin/Rscript" -e "cat(RcppParallel::CxxFlags())"`
 RCPP_PARALLEL_LIBS=`"${R_HOME}/bin/Rscript" -e "cat(RcppParallel::LdFlags())"`
 
+# Generate sources.mk using R
+"${R_HOME}/bin/Rscript" R/generate_makevars_sources.R > src/sources.mk
+
 # Substitute into Makevars
 sed -e "s|@RCPP_PARALLEL_CPPFLAGS@|${RCPP_PARALLEL_CPPFLAGS}|" \
     -e "s|@RCPP_PARALLEL_LIBS@|${RCPP_PARALLEL_LIBS}|" \
diff --git a/configure.win b/configure.win
index 080f1688..01a338bf 100755
--- a/configure.win
+++ b/configure.win
@@ -1,8 +1,12 @@
 #!/bin/sh
 
+# RcppParallel flags
 RCPP_PARALLEL_CPPFLAGS=`"${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "cat(RcppParallel::CxxFlags())"`
 RCPP_PARALLEL_LIBS=`"${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "cat(RcppParallel::LdFlags())"`
 
+# Generate sources.mk using R
+"${R_HOME}/bin/Rscript" R/generate_makevars_sources.R > src/sources.mk
+
 # Substitute into Makevars.win
 sed -e "s|@RCPP_PARALLEL_CPPFLAGS@|${RCPP_PARALLEL_CPPFLAGS}|" \
     -e "s|@RCPP_PARALLEL_LIBS@|${RCPP_PARALLEL_LIBS}|" \
diff --git a/src/Makevars.in b/src/Makevars.in
index 6a26be9a..ea013a0f 100644
--- a/src/Makevars.in
+++ b/src/Makevars.in
@@ -1,5 +1,9 @@
 CXX_STD = CXX20
 
-PKG_CPPFLAGS = @RCPP_PARALLEL_CPPFLAGS@ -DARMA_NO_DEBUG
+include sources.mk
+
+OBJECTS = $(SOURCES:.cpp=.o)
+
+PKG_CPPFLAGS = @RCPP_PARALLEL_CPPFLAGS@ -DARMA_NO_DEBUG -I.
 
 PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) @RCPP_PARALLEL_LIBS@
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 2df40048..624b82d5 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -101,38 +101,6 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// get_explog_switch
-Rcpp::String get_explog_switch();
-RcppExport SEXP _bgms_get_explog_switch() {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    rcpp_result_gen = Rcpp::wrap(get_explog_switch());
-    return rcpp_result_gen;
-END_RCPP
-}
-// rcpp_ieee754_exp
-Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x);
-RcppExport SEXP _bgms_rcpp_ieee754_exp(SEXP xSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP);
-    rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_exp(x));
-    return rcpp_result_gen;
-END_RCPP
-}
-// rcpp_ieee754_log
-Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x);
-RcppExport SEXP _bgms_rcpp_ieee754_log(SEXP xSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP);
-    rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_log(x));
-    return rcpp_result_gen;
-END_RCPP
-}
 // sample_omrf_gibbs
 IntegerMatrix sample_omrf_gibbs(int no_states, int no_variables, IntegerVector no_categories, NumericMatrix interactions, NumericMatrix thresholds, int iter);
 RcppExport SEXP _bgms_sample_omrf_gibbs(SEXP no_statesSEXP, SEXP no_variablesSEXP, SEXP no_categoriesSEXP, SEXP interactionsSEXP, SEXP thresholdsSEXP, SEXP iterSEXP) {
@@ -185,9 +153,6 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_bgms_run_bgmCompare_parallel", (DL_FUNC) &_bgms_run_bgmCompare_parallel, 36},
     {"_bgms_run_bgm_parallel", (DL_FUNC) &_bgms_run_bgm_parallel, 34},
-    {"_bgms_get_explog_switch", (DL_FUNC) &_bgms_get_explog_switch, 0},
-    {"_bgms_rcpp_ieee754_exp", (DL_FUNC) &_bgms_rcpp_ieee754_exp, 1},
-    {"_bgms_rcpp_ieee754_log", (DL_FUNC) &_bgms_rcpp_ieee754_log, 1},
     {"_bgms_sample_omrf_gibbs", (DL_FUNC) &_bgms_sample_omrf_gibbs, 6},
     {"_bgms_sample_bcomrf_gibbs", (DL_FUNC) &_bgms_sample_bcomrf_gibbs, 8},
     {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4},
diff --git a/src/bgm_helper.cpp b/src/bgm/bgm_helper.cpp
similarity index 99%
rename from src/bgm_helper.cpp
rename to src/bgm/bgm_helper.cpp
index 8f37cda9..a0696518 100644
--- a/src/bgm_helper.cpp
+++ b/src/bgm/bgm_helper.cpp
@@ -1,6 +1,6 @@
 #include <RcppArmadillo.h>
-#include "bgm_helper.h"
-#include "common_helpers.h"
+#include "bgm/bgm_helper.h"
+#include "utils/common_helpers.h"
 
 
 
diff --git a/src/bgm_helper.h b/src/bgm/bgm_helper.h
similarity index 98%
rename from src/bgm_helper.h
rename to src/bgm/bgm_helper.h
index 0cafba85..fe6e890b 100644
--- a/src/bgm_helper.h
+++ b/src/bgm/bgm_helper.h
@@ -1,7 +1,7 @@
 #pragma once
 
 #include <RcppArmadillo.h>
-#include "rng_utils.h"
+#include "rng/rng_utils.h"
 
 // Vectorize main_effect matrix
 arma::vec vectorize_main_effects_bgm(
diff --git a/src/bgm_logp_and_grad.cpp b/src/bgm/bgm_logp_and_grad.cpp
similarity index 99%
rename from src/bgm_logp_and_grad.cpp
rename to src/bgm/bgm_logp_and_grad.cpp
index a574b75e..c0ef86da 100644
--- a/src/bgm_logp_and_grad.cpp
+++ b/src/bgm/bgm_logp_and_grad.cpp
@@ -1,9 +1,9 @@
 #include <RcppArmadillo.h>
-#include "bgm_helper.h"
-#include "bgm_logp_and_grad.h"
-#include "common_helpers.h"
-#include "explog_switch.h"
-#include "variable_helpers.h"
+#include "bgm/bgm_helper.h"
+#include "bgm/bgm_logp_and_grad.h"
+#include "utils/common_helpers.h"
+#include "math/explog_switch.h"
+#include "utils/variable_helpers.h"
 
 
 
diff --git a/src/bgm_logp_and_grad.h b/src/bgm/bgm_logp_and_grad.h
similarity index 100%
rename from src/bgm_logp_and_grad.h
rename to src/bgm/bgm_logp_and_grad.h
diff --git a/src/bgm/bgm_output.h b/src/bgm/bgm_output.h
new file mode 100644
index 00000000..0768945c
--- /dev/null
+++ b/src/bgm/bgm_output.h
@@ -0,0 +1,21 @@
+#pragma once
+#include <RcppArmadillo.h>
+
+struct bgmOutput {
+  // required
+  arma::mat main_samples;
+  arma::mat pairwise_samples;
+
+  // optional (only if edge_selection)
+  arma::imat indicator_samples;
+  arma::imat allocation_samples;
+
+  // optional (only if NUTS)
+  arma::ivec treedepth_samples;
+  arma::ivec divergent_samples;
+  arma::vec energy_samples;
+
+  // metadata
+  int  chain_id = -1;
+  bool userInterrupt = false;
+};
diff --git a/src/bgm_sampler.cpp b/src/bgm/bgm_sampler.cpp
similarity index 98%
rename from src/bgm_sampler.cpp
rename to src/bgm/bgm_sampler.cpp
index 683bbf3b..680386bc 100644
--- a/src/bgm_sampler.cpp
+++ b/src/bgm/bgm_sampler.cpp
@@ -1,18 +1,19 @@
 #include <RcppArmadillo.h>
-#include "bgm_helper.h"
-#include "bgm_logp_and_grad.h"
-#include "bgm_sampler.h"
-#include "common_helpers.h"
-#include "mcmc_adaptation.h"
-#include "mcmc_hmc.h"
-#include "mcmc_leapfrog.h"
-#include "mcmc_nuts.h"
-#include "mcmc_rwm.h"
-#include "mcmc_utils.h"
-#include "sbm_edge_prior.h"
-#include "rng_utils.h"
-#include "progress_manager.h"
-#include "chainResults.h"
+#include "bgm/bgm_helper.h"
+#include "bgm/bgm_logp_and_grad.h"
+#include "bgm/bgm_sampler.h"
+#include "bgm/bgm_output.h"
+#include "mcmc/mcmc_adaptation.h"
+#include "mcmc/mcmc_hmc.h"
+#include "mcmc/mcmc_leapfrog.h"
+#include "mcmc/mcmc_nuts.h"
+#include "mcmc/mcmc_rwm.h"
+#include "mcmc/mcmc_utils.h"
+#include "priors/sbm_edge_prior.h"
+#include "sbm_edge_prior_interface.h"
+#include "rng/rng_utils.h"
+#include "utils/common_helpers.h"
+#include "utils/progress_manager.h"
 
 
 
@@ -1186,8 +1187,8 @@ void gibbs_update_step_bgm (
  *  - Parallel execution across chains is handled by `run_bgm_parallel()`;
  *    this function is for one chain only.
  */
-void run_gibbs_sampler_bgm(
-    ChainResult& chain_result,
+bgmOutput run_gibbs_sampler_bgm(
+    int chain_id,
     arma::imat observations,
     const arma::ivec& num_categories,
     const double pairwise_scale,
@@ -1221,9 +1222,6 @@ void run_gibbs_sampler_bgm(
     SafeRNG& rng,
     ProgressManager& pm
 ) {
-
-  int chain_id = chain_result.chain_id;
-
   // --- Setup: dimensions and storage structures
   const int num_variables = observations.n_cols;
   const int num_persons = observations.n_rows;
@@ -1432,22 +1430,22 @@ void run_gibbs_sampler_bgm(
     }
   }
 
+  bgmOutput chain_result;
+  chain_result.chain_id = chain_id;
   chain_result.userInterrupt = userInterrupt;
 
-  chain_result.main_effect_samples = main_effect_samples;
-  chain_result.pairwise_effect_samples = pairwise_effect_samples;
-
+  chain_result.main_samples = main_effect_samples;
+  chain_result.pairwise_samples = pairwise_effect_samples;
   if (update_method == nuts) {
     chain_result.treedepth_samples = treedepth_samples;
     chain_result.divergent_samples = divergent_samples;
     chain_result.energy_samples    = energy_samples;
   }
-
   if (edge_selection) {
     chain_result.indicator_samples = indicator_samples;
-
     if (edge_prior == Stochastic_Block)
       chain_result.allocation_samples = allocation_samples;
   }
 
+  return chain_result;
 }
\ No newline at end of file
diff --git a/src/bgm_sampler.h b/src/bgm/bgm_sampler.h
similarity index 91%
rename from src/bgm_sampler.h
rename to src/bgm/bgm_sampler.h
index 2fcbf357..3e280ee1 100644
--- a/src/bgm_sampler.h
+++ b/src/bgm/bgm_sampler.h
@@ -1,13 +1,14 @@
 #pragma once
 #include <RcppArmadillo.h>
-#include "common_helpers.h"
+#include "utils/common_helpers.h"
+#include "bgm/bgm_output.h"
+
 // forward declaration
 struct SafeRNG;
 class ProgressManager;
-struct ChainResult;
 
-void run_gibbs_sampler_bgm(
-    ChainResult& chain_result,
+bgmOutput run_gibbs_sampler_bgm(
+    int chain_id,
     arma::imat observations,
     const arma::ivec& num_categories,
     const double pairwise_scale,
diff --git a/src/bgmCompare_helper.cpp b/src/bgmCompare/bgmCompare_helper.cpp
similarity index 99%
rename from src/bgmCompare_helper.cpp
rename to src/bgmCompare/bgmCompare_helper.cpp
index f73fc5fb..10706afa 100644
--- a/src/bgmCompare_helper.cpp
+++ b/src/bgmCompare/bgmCompare_helper.cpp
@@ -1,7 +1,7 @@
 #include <RcppArmadillo.h>
 #include <cmath>
-#include "bgmCompare_helper.h"
-#include "common_helpers.h"
+#include "bgmCompare/bgmCompare_helper.h"
+#include "utils/common_helpers.h"
 
 
 
diff --git a/src/bgmCompare_helper.h b/src/bgmCompare/bgmCompare_helper.h
similarity index 99%
rename from src/bgmCompare_helper.h
rename to src/bgmCompare/bgmCompare_helper.h
index 491609fd..05759b3c 100644
--- a/src/bgmCompare_helper.h
+++ b/src/bgmCompare/bgmCompare_helper.h
@@ -1,7 +1,7 @@
 #pragma once
 
 #include <RcppArmadillo.h>
-#include "rng_utils.h"
+#include "rng/rng_utils.h"
 
 
 
diff --git a/src/bgmCompare_logp_and_grad.cpp b/src/bgmCompare/bgmCompare_logp_and_grad.cpp
similarity index 99%
rename from src/bgmCompare_logp_and_grad.cpp
rename to src/bgmCompare/bgmCompare_logp_and_grad.cpp
index 4807ff62..5b4c15b6 100644
--- a/src/bgmCompare_logp_and_grad.cpp
+++ b/src/bgmCompare/bgmCompare_logp_and_grad.cpp
@@ -1,10 +1,10 @@
 #include <RcppArmadillo.h>
-#include "bgmCompare_helper.h"
-#include "bgmCompare_logp_and_grad.h"
+#include "bgmCompare/bgmCompare_helper.h"
+#include "bgmCompare/bgmCompare_logp_and_grad.h"
 #include <cmath>
-#include "explog_switch.h"
-#include "common_helpers.h"
-#include "variable_helpers.h"
+#include "math/explog_switch.h"
+#include "utils/common_helpers.h"
+#include "utils/variable_helpers.h"
 
 
 
diff --git a/src/bgmCompare_logp_and_grad.h b/src/bgmCompare/bgmCompare_logp_and_grad.h
similarity index 100%
rename from src/bgmCompare_logp_and_grad.h
rename to src/bgmCompare/bgmCompare_logp_and_grad.h
diff --git a/src/sampler_output.h b/src/bgmCompare/bgmCompare_output.h
similarity index 97%
rename from src/sampler_output.h
rename to src/bgmCompare/bgmCompare_output.h
index 8f899d00..148a1984 100644
--- a/src/sampler_output.h
+++ b/src/bgmCompare/bgmCompare_output.h
@@ -21,13 +21,15 @@
  *  - chain_id:         Identifier of the chain.
  *  - has_indicator:    True if indicator samples are stored.
  */
-struct SamplerOutput {
+struct bgmCompareOutput {
   arma::mat main_samples;
   arma::mat pairwise_samples;
   arma::imat indicator_samples;
+
   arma::ivec treedepth_samples;
   arma::ivec divergent_samples;
   arma::vec energy_samples;
+
   int chain_id;
   bool has_indicator;
   bool userInterrupt;
diff --git a/src/bgmCompare_sampler.cpp b/src/bgmCompare/bgmCompare_sampler.cpp
similarity index 99%
rename from src/bgmCompare_sampler.cpp
rename to src/bgmCompare/bgmCompare_sampler.cpp
index 14853465..cd955e32 100644
--- a/src/bgmCompare_sampler.cpp
+++ b/src/bgmCompare/bgmCompare_sampler.cpp
@@ -1,19 +1,19 @@
 #include <RcppArmadillo.h>
-#include "bgmCompare_helper.h"
-#include "bgmCompare_logp_and_grad.h"
-#include "bgmCompare_sampler.h"
-#include "common_helpers.h"
-#include "mcmc_adaptation.h"
-#include "mcmc_hmc.h"
-#include "mcmc_leapfrog.h"
-#include "mcmc_nuts.h"
-#include "mcmc_rwm.h"
-#include "mcmc_utils.h"
-#include "rng_utils.h"
-#include "sampler_output.h"
-#include "explog_switch.h"
+#include "bgmCompare/bgmCompare_helper.h"
+#include "bgmCompare/bgmCompare_logp_and_grad.h"
+#include "bgmCompare/bgmCompare_sampler.h"
+#include "bgmCompare/bgmCompare_output.h"
+#include "mcmc/mcmc_adaptation.h"
+#include "mcmc/mcmc_hmc.h"
+#include "mcmc/mcmc_leapfrog.h"
+#include "mcmc/mcmc_nuts.h"
+#include "mcmc/mcmc_rwm.h"
+#include "mcmc/mcmc_utils.h"
+#include "rng/rng_utils.h"
+#include "math/explog_switch.h"
 #include <string>
-#include "progress_manager.h"
+#include "utils/progress_manager.h"
+#include "utils/common_helpers.h"
 
 
 
@@ -1564,7 +1564,7 @@ void gibbs_update_step_bgmcompare (
  *  - This function runs entirely in C++ and is wrapped for parallel execution
  *    via `GibbsCompareChainRunner`.
  */
-SamplerOutput run_gibbs_sampler_bgmCompare(
+bgmCompareOutput run_gibbs_sampler_bgmCompare(
     int chain_id,
     arma::imat observations,
     const int num_groups,
@@ -1761,7 +1761,7 @@ SamplerOutput run_gibbs_sampler_bgmCompare(
     }
   }
 
-  SamplerOutput out;
+  bgmCompareOutput out;
   out.chain_id = chain_id;
   out.main_samples = main_effect_samples;
   out.pairwise_samples = pairwise_effect_samples;
@@ -1774,7 +1774,7 @@ SamplerOutput run_gibbs_sampler_bgmCompare(
   } else {
     out.indicator_samples = arma::imat();
   }
- out.userInterrupt = userInterrupt;
+  out.userInterrupt = userInterrupt;
 
   return out;
 }
\ No newline at end of file
diff --git a/src/bgmCompare_sampler.h b/src/bgmCompare/bgmCompare_sampler.h
similarity index 91%
rename from src/bgmCompare_sampler.h
rename to src/bgmCompare/bgmCompare_sampler.h
index 781d7f0d..83bcdbaa 100644
--- a/src/bgmCompare_sampler.h
+++ b/src/bgmCompare/bgmCompare_sampler.h
@@ -1,14 +1,14 @@
 #pragma once
 
 #include <RcppArmadillo.h>
-#include "common_helpers.h"
+#include "utils/common_helpers.h"
+#include "bgmCompare/bgmCompare_output.h"
 #include <string>
 
-struct SamplerOutput;
 struct SafeRNG;
 class ProgressManager;
 
-SamplerOutput run_gibbs_sampler_bgmCompare(
+bgmCompareOutput run_gibbs_sampler_bgmCompare(
     int chain_id,
     arma::imat observations,
     const int num_groups,
diff --git a/src/bgmCompare_parallel.cpp b/src/bgmCompare_interface.cpp
similarity index 95%
rename from src/bgmCompare_parallel.cpp
rename to src/bgmCompare_interface.cpp
index 24464023..5b0fa983 100644
--- a/src/bgmCompare_parallel.cpp
+++ b/src/bgmCompare_interface.cpp
@@ -1,15 +1,15 @@
 // [[Rcpp::depends(RcppParallel, RcppArmadillo, dqrng)]]
 #include <RcppArmadillo.h>
-#include "bgmCompare_sampler.h"
-#include "rng_utils.h" // must be included before RcppParallel
+#include "bgmCompare/bgmCompare_sampler.h"
+#include "rng/rng_utils.h" // must be included before RcppParallel
 #include <RcppParallel.h>
 #include <tbb/global_control.h>
 #include <vector>
 #include <string>
-#include "progress_manager.h"
-#include "sampler_output.h"
-#include "mcmc_adaptation.h"
-#include "common_helpers.h"
+#include "utils/progress_manager.h"
+#include "bgmCompare/bgmCompare_output.h"
+#include "mcmc/mcmc_adaptation.h"
+#include "utils/common_helpers.h"
 
 using namespace RcppParallel;
 
@@ -22,24 +22,19 @@ using namespace RcppParallel;
  *  - error: True if the chain terminated with an error, false otherwise.
  *  - error_msg: Error message if an error occurred (empty if none).
  *  - chain_id: Integer identifier for the chain (1-based).
- *  - result: SamplerOutput object containing chain results
+ *  - result: bgmCompareOutput object containing chain results
  *    (samples, diagnostics, metadata).
  *
  * Usage:
  *  - Used in parallel execution to collect results from each chain.
  *  - Checked after execution to propagate errors or assemble outputs
  *    into an R-accessible list.
- *
- * Notes:
- *  - This struct mirrors `ChainResult` from bgm, but stores a
- *    `SamplerOutput` instead of an `Rcpp::List`.
- *  - If `error == true`, the `result` field should be ignored.
  */
-struct ChainResultCompare {
+struct bgmCompareChainResult {
   bool error;
   std::string error_msg;
   int chain_id;
-  SamplerOutput result;
+  bgmCompareOutput result;
 };
 
 
@@ -90,7 +85,7 @@ struct ChainResultCompare {
  *  - hmc_num_leapfrogs: Number of leapfrog steps (HMC).
  *
  * Output:
- *  - results: Vector of `ChainResultCompare` objects, one per chain, filled in place.
+ *  - results: Vector of `bgmCompareChainResult` objects, one per chain, filled in place.
  *
  * Notes:
  *  - Each worker instance is shared across threads but invoked with different
@@ -135,7 +130,7 @@ struct GibbsCompareChainRunner : public Worker {
   const int hmc_num_leapfrogs;
   ProgressManager& pm;
   // output
-  std::vector<ChainResultCompare>& results;
+  std::vector<bgmCompareChainResult>& results;
 
   GibbsCompareChainRunner(
     const arma::imat& observations_master,
@@ -172,7 +167,7 @@ struct GibbsCompareChainRunner : public Worker {
     const UpdateMethod update_method,
     const int hmc_num_leapfrogs,
     ProgressManager& pm,
-    std::vector<ChainResultCompare>& results
+    std::vector<bgmCompareChainResult>& results
   ) :
     observations_master(observations_master),
     num_groups(num_groups),
@@ -213,7 +208,7 @@ struct GibbsCompareChainRunner : public Worker {
 
   void operator()(std::size_t begin, std::size_t end) {
     for (std::size_t i = begin; i < end; ++i) {
-      ChainResultCompare out;
+      bgmCompareChainResult out;
       out.chain_id = static_cast<int>(i + 1);
       out.error = false;
 
@@ -229,7 +224,7 @@ struct GibbsCompareChainRunner : public Worker {
         arma::imat observations = observations_master;
 
         // run sampler (pure C++)
-        SamplerOutput result = run_gibbs_sampler_bgmCompare(
+        bgmCompareOutput result = run_gibbs_sampler_bgmCompare(
           out.chain_id,
           observations,
           num_groups,
@@ -386,7 +381,7 @@ Rcpp::List run_bgmCompare_parallel(
     int hmc_num_leapfrogs,
     int progress_type
 ) {
-  std::vector<ChainResultCompare> results(num_chains);
+  std::vector<bgmCompareChainResult> results(num_chains);
 
   // per-chain seeds
   std::vector<SafeRNG> chain_rngs(num_chains);
diff --git a/src/bgm_parallel.cpp b/src/bgm_interface.cpp
similarity index 84%
rename from src/bgm_parallel.cpp
rename to src/bgm_interface.cpp
index f0e6c941..86834b16 100644
--- a/src/bgm_parallel.cpp
+++ b/src/bgm_interface.cpp
@@ -1,20 +1,43 @@
 // [[Rcpp::depends(RcppParallel, RcppArmadillo, dqrng)]]
 #include <RcppArmadillo.h>
-#include "rng_utils.h" // must be included before RcppParallel
+#include "rng/rng_utils.h" // must be included before RcppParallel
 #include <RcppParallel.h>
-#include "bgm_sampler.h"
+#include "bgm/bgm_sampler.h"
+#include "bgm/bgm_output.h"
 #include <tbb/global_control.h>
 #include <vector>
 #include <string>
-#include "progress_manager.h"
-#include "mcmc_adaptation.h"
-#include "common_helpers.h"
-#include "chainResults.h"
+#include "mcmc/mcmc_adaptation.h"
+#include "utils/progress_manager.h"
+#include "utils/common_helpers.h"
 
 using namespace RcppParallel;
 
 
 
+/**
+ * Container for the result of a single MCMC chain (bgm model).
+ *
+ * Fields:
+ *  - error: True if the chain terminated with an error, false otherwise.
+ *  - error_msg: Error message if an error occurred (empty if none).
+ *  - chain_id: Integer identifier for the chain (1-based).
+ *  - result: bgmOutput object containing chain results
+ *    (samples, diagnostics, metadata).
+ *
+ * Usage:
+ *  - Used in parallel execution to collect results from each chain.
+ *  - Checked after execution to propagate errors or assemble outputs
+ *    into an R-accessible list.
+ */
+struct bgmChainResult {
+  bool error;
+  std::string error_msg;
+  int chain_id;
+  bgmOutput result;
+};
+
+
 
 /**
  * Worker struct for running a single Gibbs sampling chain in parallel (bgm model).
@@ -40,7 +63,7 @@ using namespace RcppParallel;
 struct GibbsChainRunner : public Worker {
   const arma::imat& observations;
   const arma::ivec& num_categories;
-  double  pairwise_scale;
+  double pairwise_scale;
   const EdgePrior edge_prior;
   const arma::mat& inclusion_probability;
   double beta_bernoulli_alpha;
@@ -74,12 +97,12 @@ struct GibbsChainRunner : public Worker {
   ProgressManager& pm;
 
   // output buffer
-  std::vector<ChainResult>& results;
+  std::vector<bgmChainResult>& results;
 
   GibbsChainRunner(
     const arma::imat& observations,
     const arma::ivec& num_categories,
-    double  pairwise_scale,
+    double pairwise_scale,
     const EdgePrior edge_prior,
     const arma::mat& inclusion_probability,
     double beta_bernoulli_alpha,
@@ -109,11 +132,11 @@ struct GibbsChainRunner : public Worker {
     bool learn_mass_matrix,
     const std::vector<SafeRNG>& chain_rngs,
     ProgressManager& pm,
-    std::vector<ChainResult>& results
+    std::vector<bgmChainResult>& results
   ) :
     observations(observations),
     num_categories(num_categories),
-     pairwise_scale( pairwise_scale),
+    pairwise_scale( pairwise_scale),
     edge_prior(edge_prior),
     inclusion_probability(inclusion_probability),
     beta_bernoulli_alpha(beta_bernoulli_alpha),
@@ -148,16 +171,16 @@ struct GibbsChainRunner : public Worker {
 
   void operator()(std::size_t begin, std::size_t end) {
     for (std::size_t i = begin; i < end; ++i) {
-
-      ChainResult& chain_result = results[i];
+      bgmChainResult chain_result;
       chain_result.chain_id = static_cast<int>(i + 1);
       chain_result.error = false;
-      SafeRNG rng = chain_rngs[i];
 
       try {
+        // per-chain RNG
+        SafeRNG rng = chain_rngs[i];
 
-        run_gibbs_sampler_bgm(
-          chain_result,
+        bgmOutput result = run_gibbs_sampler_bgm(
+          chain_result.chain_id,
           observations,
           num_categories,
           pairwise_scale,
@@ -192,6 +215,8 @@ struct GibbsChainRunner : public Worker {
           pm
         );
 
+        chain_result.result = result;
+
       } catch (std::exception& e) {
         chain_result.error = true;
         chain_result.error_msg = e.what();
@@ -199,6 +224,8 @@ struct GibbsChainRunner : public Worker {
         chain_result.error = true;
         chain_result.error_msg = "Unknown error";
       }
+
+      results[i] = chain_result;
     }
   }
 };
@@ -288,7 +315,7 @@ Rcpp::List run_bgm_parallel(
     int seed,
     int progress_type
 ) {
-  std::vector<ChainResult> results(num_chains);
+  std::vector<bgmChainResult> results(num_chains);
 
   // Prepare one independent RNG per chain via jump()
   std::vector<SafeRNG> chain_rngs(num_chains);
@@ -328,31 +355,32 @@ Rcpp::List run_bgm_parallel(
         Rcpp::Named("chain_id") = results[i].chain_id
       );
     } else {
-        Rcpp::List chain_i;
-        chain_i["main_samples"] = results[i].main_effect_samples;
-        chain_i["pairwise_samples"] = results[i].pairwise_effect_samples;
+      const auto& r = results[i].result;
+      Rcpp::List chain_i;
+      chain_i["main_samples"] = r.main_samples;
+      chain_i["pairwise_samples"] = r.pairwise_samples;
 
-        if (update_method_enum == nuts) {
-          chain_i["treedepth__"] = results[i].treedepth_samples;
-          chain_i["divergent__"] = results[i].divergent_samples;
-          chain_i["energy__"] = results[i].energy_samples;
-        }
+      if (update_method_enum == nuts) {
+        chain_i["treedepth__"] = r.treedepth_samples;
+        chain_i["divergent__"] = r.divergent_samples;
+        chain_i["energy__"] = r.energy_samples;
+      }
 
-        if (edge_selection) {
-          chain_i["indicator_samples"] = results[i].indicator_samples;
+      if (edge_selection) {
+        chain_i["indicator_samples"] = r.indicator_samples;
 
-          if (edge_prior_enum == Stochastic_Block)
-            chain_i["allocations"] = results[i].allocation_samples;
-        }
+        if (edge_prior_enum == Stochastic_Block)
+          chain_i["allocations"] = r.allocation_samples;
+      }
 
-        chain_i["userInterrupt"] = results[i].userInterrupt;
-        chain_i["chain_id"] = results[i].chain_id;
+      chain_i["userInterrupt"] = r.userInterrupt;
+      chain_i["chain_id"] = r.chain_id;
 
-        output[i] = chain_i;
+      output[i] = chain_i;
     }
   }
 
   pm.finish();
 
   return output;
-}
+}
\ No newline at end of file
diff --git a/src/chainResults.h b/src/chainResults.h
deleted file mode 100644
index 395e9d42..00000000
--- a/src/chainResults.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#pragma once
-
-#include <RcppArmadillo.h>
-#include <string>
-
-/**
- * Container for the result of a single MCMC chain.
- *
- * Fields:
- *  - error: True if the chain terminated with an error, false otherwise.
- *  - error_msg: Error message in case of failure (empty if no error).
- *  - chain_id: Integer identifier for the chain (1-based).
- *  - result: Rcpp::List containing the chain’s outputs (samples, diagnostics, etc.).
- *
- * Usage:
- *  - Used in parallel samplers to collect per-chain results.
- *  - Checked after execution to propagate errors or assemble outputs into R.
- */
-struct ChainResult {
-  bool error;
-  std::string error_msg;
-  int chain_id;
-  bool userInterrupt;
-  arma::mat main_effect_samples;
-  arma::mat pairwise_effect_samples;
-  arma::ivec treedepth_samples;
-  arma::ivec divergent_samples;
-  arma::vec energy_samples;
-  arma::imat indicator_samples;
-  arma::imat allocation_samples;
-};
diff --git a/src/data_simulation.cpp b/src/data_simulation.cpp
index 55639a2c..b4bf38c3 100644
--- a/src/data_simulation.cpp
+++ b/src/data_simulation.cpp
@@ -1,4 +1,4 @@
-#include "explog_switch.h"
+#include "math/explog_switch.h"
 #include <Rcpp.h>
 using namespace Rcpp;
 
diff --git a/src/custom_exp.cpp b/src/math/custom_exp.cpp
similarity index 85%
rename from src/custom_exp.cpp
rename to src/math/custom_exp.cpp
index e6fe3164..2abe8cb3 100644
--- a/src/custom_exp.cpp
+++ b/src/math/custom_exp.cpp
@@ -1,7 +1,6 @@
-#include "explog_switch.h"
+#include "math/explog_switch.h"
 #include "Rcpp.h"
 
-// [[Rcpp::export]]
 Rcpp::String get_explog_switch() {
 #if USE_CUSTOM_LOG
       return "custom";
@@ -10,7 +9,6 @@ Rcpp::String get_explog_switch() {
   #endif
 }
 
-// [[Rcpp::export]]
 Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x) {
   Rcpp::NumericVector y(x.size());
   for (int i = 0; i < x.size(); i++) {
@@ -19,7 +17,6 @@ Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x) {
   return y;
 }
 
-// [[Rcpp::export]]
 Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x) {
   Rcpp::NumericVector y(x.size());
   for (int i = 0; i < x.size(); i++) {
diff --git a/src/e_arma_exp.h b/src/math/e_arma_exp.h
similarity index 100%
rename from src/e_arma_exp.h
rename to src/math/e_arma_exp.h
diff --git a/src/e_exp.cpp b/src/math/e_exp.cpp
similarity index 100%
rename from src/e_exp.cpp
rename to src/math/e_exp.cpp
diff --git a/src/e_exp.h b/src/math/e_exp.h
similarity index 100%
rename from src/e_exp.h
rename to src/math/e_exp.h
diff --git a/src/explog_switch.h b/src/math/explog_switch.h
similarity index 94%
rename from src/explog_switch.h
rename to src/math/explog_switch.h
index 9c174cc4..3fce818f 100644
--- a/src/explog_switch.h
+++ b/src/math/explog_switch.h
@@ -20,8 +20,8 @@
 
 #if USE_CUSTOM_LOG
 
-#include "e_exp.h"
-#include "e_arma_exp.h"
+#include "math/e_exp.h"
+#include "math/e_arma_exp.h"
 
 #define MY_EXP __ieee754_exp
 #define MY_LOG __ieee754_log
diff --git a/src/mcmc_adaptation.h b/src/mcmc/mcmc_adaptation.h
similarity index 99%
rename from src/mcmc_adaptation.h
rename to src/mcmc/mcmc_adaptation.h
index 4aa172c9..cf6184ff 100644
--- a/src/mcmc_adaptation.h
+++ b/src/mcmc/mcmc_adaptation.h
@@ -5,9 +5,9 @@
 #include <cmath>
 #include <algorithm>
 #include <vector>
-#include "mcmc_utils.h"
-#include "mcmc_rwm.h"
-#include "explog_switch.h"
+#include "mcmc/mcmc_utils.h"
+#include "mcmc/mcmc_rwm.h"
+#include "math/explog_switch.h"
 
 class DualAveraging {
 public:
diff --git a/src/mcmc_hmc.cpp b/src/mcmc/mcmc_hmc.cpp
similarity index 89%
rename from src/mcmc_hmc.cpp
rename to src/mcmc/mcmc_hmc.cpp
index e851a24a..adb3c022 100644
--- a/src/mcmc_hmc.cpp
+++ b/src/mcmc/mcmc_hmc.cpp
@@ -1,9 +1,9 @@
 #include <RcppArmadillo.h>
 #include <functional>
-#include "mcmc_hmc.h"
-#include "mcmc_leapfrog.h"
-#include "mcmc_utils.h"
-#include "rng_utils.h"
+#include "mcmc/mcmc_hmc.h"
+#include "mcmc/mcmc_leapfrog.h"
+#include "mcmc/mcmc_utils.h"
+#include "rng/rng_utils.h"
 
 
 
diff --git a/src/mcmc_hmc.h b/src/mcmc/mcmc_hmc.h
similarity index 92%
rename from src/mcmc_hmc.h
rename to src/mcmc/mcmc_hmc.h
index 6492cd56..32e9d4a1 100644
--- a/src/mcmc_hmc.h
+++ b/src/mcmc/mcmc_hmc.h
@@ -2,7 +2,7 @@
 
 #include <RcppArmadillo.h>
 #include <functional>
-#include "mcmc_utils.h"
+#include "mcmc/mcmc_utils.h"
 struct SafeRNG;
 
 SamplerResult hmc_sampler(
diff --git a/src/mcmc_leapfrog.cpp b/src/mcmc/mcmc_leapfrog.cpp
similarity index 97%
rename from src/mcmc_leapfrog.cpp
rename to src/mcmc/mcmc_leapfrog.cpp
index 49c7d10e..ceadcfe8 100644
--- a/src/mcmc_leapfrog.cpp
+++ b/src/mcmc/mcmc_leapfrog.cpp
@@ -1,7 +1,7 @@
 #include <RcppArmadillo.h>
 #include <functional>
-#include "mcmc_leapfrog.h"
-#include "mcmc_memoization.h"
+#include "mcmc/mcmc_leapfrog.h"
+#include "mcmc/mcmc_memoization.h"
 
 
 
diff --git a/src/mcmc_leapfrog.h b/src/mcmc/mcmc_leapfrog.h
similarity index 95%
rename from src/mcmc_leapfrog.h
rename to src/mcmc/mcmc_leapfrog.h
index b69164f1..bd7a5397 100644
--- a/src/mcmc_leapfrog.h
+++ b/src/mcmc/mcmc_leapfrog.h
@@ -3,7 +3,7 @@
 
 #include <RcppArmadillo.h>
 #include <functional>
-#include "mcmc_memoization.h"
+#include "mcmc/mcmc_memoization.h"
 
 
 
diff --git a/src/mcmc_memoization.h b/src/mcmc/mcmc_memoization.h
similarity index 100%
rename from src/mcmc_memoization.h
rename to src/mcmc/mcmc_memoization.h
diff --git a/src/mcmc_nuts.cpp b/src/mcmc/mcmc_nuts.cpp
similarity index 98%
rename from src/mcmc_nuts.cpp
rename to src/mcmc/mcmc_nuts.cpp
index 96884dc6..2335fabc 100644
--- a/src/mcmc_nuts.cpp
+++ b/src/mcmc/mcmc_nuts.cpp
@@ -1,10 +1,10 @@
 #include <RcppArmadillo.h>
 #include <functional>
-#include "mcmc_leapfrog.h"
-#include "mcmc_memoization.h"
-#include "mcmc_nuts.h"
-#include "mcmc_utils.h"
-#include "rng_utils.h"
+#include "mcmc/mcmc_leapfrog.h"
+#include "mcmc/mcmc_memoization.h"
+#include "mcmc/mcmc_nuts.h"
+#include "mcmc/mcmc_utils.h"
+#include "rng/rng_utils.h"
 
 
 /**
diff --git a/src/mcmc_nuts.h b/src/mcmc/mcmc_nuts.h
similarity index 94%
rename from src/mcmc_nuts.h
rename to src/mcmc/mcmc_nuts.h
index 613ffc2b..f678cf4b 100644
--- a/src/mcmc_nuts.h
+++ b/src/mcmc/mcmc_nuts.h
@@ -3,9 +3,9 @@
 #include <RcppArmadillo.h>
 #include <functional>
 #include <utility>
-#include "mcmc_leapfrog.h"
-#include "mcmc_memoization.h"
-#include "mcmc_utils.h"
+#include "mcmc/mcmc_leapfrog.h"
+#include "mcmc/mcmc_memoization.h"
+#include "mcmc/mcmc_utils.h"
 struct SafeRNG;
 
 
diff --git a/src/mcmc_rwm.cpp b/src/mcmc/mcmc_rwm.cpp
similarity index 93%
rename from src/mcmc_rwm.cpp
rename to src/mcmc/mcmc_rwm.cpp
index a91dc8f4..a07ad105 100644
--- a/src/mcmc_rwm.cpp
+++ b/src/mcmc/mcmc_rwm.cpp
@@ -1,9 +1,9 @@
 #include <RcppArmadillo.h>
 #include <cmath>
 #include <functional>
-#include "mcmc_utils.h"
-#include "mcmc_rwm.h"
-#include "rng_utils.h"
+#include "mcmc/mcmc_utils.h"
+#include "mcmc/mcmc_rwm.h"
+#include "rng/rng_utils.h"
 
 
 /**
diff --git a/src/mcmc_rwm.h b/src/mcmc/mcmc_rwm.h
similarity index 89%
rename from src/mcmc_rwm.h
rename to src/mcmc/mcmc_rwm.h
index 73227642..fcce0254 100644
--- a/src/mcmc_rwm.h
+++ b/src/mcmc/mcmc_rwm.h
@@ -3,7 +3,7 @@
 #include <RcppArmadillo.h>
 #include <cmath>
 #include <functional>
-#include "mcmc_utils.h"
+#include "mcmc/mcmc_utils.h"
 struct SafeRNG;
 
 
diff --git a/src/mcmc_utils.cpp b/src/mcmc/mcmc_utils.cpp
similarity index 96%
rename from src/mcmc_utils.cpp
rename to src/mcmc/mcmc_utils.cpp
index 4c2328c4..25c616c4 100644
--- a/src/mcmc_utils.cpp
+++ b/src/mcmc/mcmc_utils.cpp
@@ -1,9 +1,9 @@
 #include <RcppArmadillo.h>
 #include <cmath>
 #include <functional>
-#include "mcmc_leapfrog.h"
-#include "mcmc_utils.h"
-#include "rng_utils.h"
+#include "mcmc/mcmc_leapfrog.h"
+#include "mcmc/mcmc_utils.h"
+#include "rng/rng_utils.h"
 
 
 
diff --git a/src/mcmc_utils.h b/src/mcmc/mcmc_utils.h
similarity index 99%
rename from src/mcmc_utils.h
rename to src/mcmc/mcmc_utils.h
index 9c935706..c84c0393 100644
--- a/src/mcmc_utils.h
+++ b/src/mcmc/mcmc_utils.h
@@ -5,7 +5,7 @@
 #include <cmath>
 #include <functional>
 #include <memory>
-#include "explog_switch.h"
+#include "math/explog_switch.h"
 struct SafeRNG;
 
 // (only if <algorithm> didn’t already provide it under C++17)
diff --git a/src/sbm_edge_prior.cpp b/src/priors/sbm_edge_prior.cpp
similarity index 91%
rename from src/sbm_edge_prior.cpp
rename to src/priors/sbm_edge_prior.cpp
index 182201bc..b24d2f57 100644
--- a/src/sbm_edge_prior.cpp
+++ b/src/priors/sbm_edge_prior.cpp
@@ -1,6 +1,6 @@
 #include <RcppArmadillo.h>
-#include "rng_utils.h"
-#include "explog_switch.h"
+#include "rng/rng_utils.h"
+#include "math/explog_switch.h"
 
 // ----------------------------------------------------------------------------|
 // The c++ code below is based on the R code accompanying the paper:
@@ -58,30 +58,6 @@ arma::mat add_row_col_block_prob_matrix(arma::mat X,
 }
 
 
-// ----------------------------------------------------------------------------|
-// Compute partition coefficient for the MFM - SBM
-// ----------------------------------------------------------------------------|
-// [[Rcpp::export]]
-arma::vec compute_Vn_mfm_sbm(arma::uword no_variables,
-                             double dirichlet_alpha,
-                             arma::uword t_max,
-                             double lambda) {
-  arma::vec log_Vn(t_max);
-  double r;
-
-  for(arma::uword t = 0; t < t_max; t++) {
-    r = -INFINITY; // initialize log-coefficient at -Inf
-    for(arma::uword k = t; k <= 500; k++){
-      arma::vec b_linspace_1 = arma::linspace(k-t+1,k+1,t+1); // numerator = b*(b-1)*...*(b-|C|+1)
-      arma::vec b_linspace_2 = arma::linspace((k+1)*dirichlet_alpha,(k+1)*dirichlet_alpha+no_variables-1, no_variables); // denominator b*e*(b*e+1)*...*(b*e+p-1)
-      double b = arma::accu(arma::log(b_linspace_1))-arma::accu(arma::log(b_linspace_2)) + R::dpois((k+1)-1, lambda, true); // sum(log(numerator)) - sum(log(denominator)) + log(P=(k+1|lambda))
-      double m = std::max(b,r);  // scaling factor for log-sum-exp formula
-      r = MY_LOG(MY_EXP(r-m) +  MY_EXP(b-m)) + m; // update r using log-sum-exp formula to ensure numerical stability and avoid underflow
-    }
-    log_Vn(t) = r;
-  }
-  return log_Vn;
-}
 
 // ----------------------------------------------------------------------------|
 // Compute log-likelihood for the MFM - SBM
diff --git a/src/sbm_edge_prior.h b/src/priors/sbm_edge_prior.h
similarity index 80%
rename from src/sbm_edge_prior.h
rename to src/priors/sbm_edge_prior.h
index a782347b..b4a1821f 100644
--- a/src/sbm_edge_prior.h
+++ b/src/priors/sbm_edge_prior.h
@@ -4,13 +4,6 @@
 struct SafeRNG;
 
 
-// ----------------------------------------------------------------------------|
-// Compute partition coefficient for the MFM - SBM
-// ----------------------------------------------------------------------------|
-arma::vec compute_Vn_mfm_sbm(arma::uword no_variables,
-                                        double dirichlet_alpha,
-                                        arma::uword t_max,
-                                        double lambda);
 
 // ----------------------------------------------------------------------------|
 // Sample the block allocations for the MFM - SBM
diff --git a/src/rng_utils.h b/src/rng/rng_utils.h
similarity index 100%
rename from src/rng_utils.h
rename to src/rng/rng_utils.h
diff --git a/src/sbm_edge_prior_interface.cpp b/src/sbm_edge_prior_interface.cpp
new file mode 100644
index 00000000..8e040b90
--- /dev/null
+++ b/src/sbm_edge_prior_interface.cpp
@@ -0,0 +1,34 @@
+#include <RcppArmadillo.h>
+#include "math/explog_switch.h"
+
+// ----------------------------------------------------------------------------|
+// The c++ code below is based on the R code accompanying the paper:
+//  Geng, J., Bhattacharya, A., & Pati, D. (2019). Probabilistic Community
+//  Detection With Unknown Number of Communities, Journal of the American
+//  Statistical Association, 114:526, 893-905, DOI:10.1080/01621459.2018.1458618
+// ----------------------------------------------------------------------------|
+
+// ----------------------------------------------------------------------------|
+// Compute partition coefficient for the MFM - SBM
+// ----------------------------------------------------------------------------|
+// [[Rcpp::export]]
+arma::vec compute_Vn_mfm_sbm(arma::uword no_variables,
+                             double dirichlet_alpha,
+                             arma::uword t_max,
+                             double lambda) {
+  arma::vec log_Vn(t_max);
+  double r;
+
+  for(arma::uword t = 0; t < t_max; t++) {
+    r = -INFINITY; // initialize log-coefficient at -Inf
+    for(arma::uword k = t; k <= 500; k++){
+      arma::vec b_linspace_1 = arma::linspace(k-t+1,k+1,t+1); // numerator = b*(b-1)*...*(b-|C|+1)
+      arma::vec b_linspace_2 = arma::linspace((k+1)*dirichlet_alpha,(k+1)*dirichlet_alpha+no_variables-1, no_variables); // denominator b*e*(b*e+1)*...*(b*e+p-1)
+      double b = arma::accu(ARMA_MY_LOG(b_linspace_1))-arma::accu(ARMA_MY_LOG(b_linspace_2)) + R::dpois((k+1)-1, lambda, true); // sum(log(numerator)) - sum(log(denominator)) + log(P=(k+1|lambda))
+      double m = std::max(b,r);  // scaling factor for log-sum-exp formula
+      r = MY_LOG(MY_EXP(r-m) +  MY_EXP(b-m)) + m; // update r using log-sum-exp formula to ensure numerical stability and avoid underflow
+    }
+    log_Vn(t) = r;
+  }
+  return log_Vn;
+}
\ No newline at end of file
diff --git a/src/sbm_edge_prior_interface.h b/src/sbm_edge_prior_interface.h
new file mode 100644
index 00000000..ac746a90
--- /dev/null
+++ b/src/sbm_edge_prior_interface.h
@@ -0,0 +1,9 @@
+#include <RcppArmadillo.h>
+#include "math/explog_switch.h"
+
+
+
+arma::vec compute_Vn_mfm_sbm(arma::uword no_variables,
+                             double dirichlet_alpha,
+                             arma::uword t_max,
+                             double lambda);
\ No newline at end of file
diff --git a/src/sources.mk b/src/sources.mk
new file mode 100644
index 00000000..8a780ef0
--- /dev/null
+++ b/src/sources.mk
@@ -0,0 +1,27 @@
+# ------------------------------------------------------------------
+# THIS FILE IS AUTO-GENERATED - DO NOT EDIT
+# Generated by configure
+# To add C++ code, place .cpp files anywhere under src/
+# ------------------------------------------------------------------
+SOURCES = \
+  RcppExports.cpp \
+  bgm/bgm_helper.cpp \
+  bgm/bgm_logp_and_grad.cpp \
+  bgm/bgm_sampler.cpp \
+  bgmCompare/bgmCompare_helper.cpp \
+  bgmCompare/bgmCompare_logp_and_grad.cpp \
+  bgmCompare/bgmCompare_sampler.cpp \
+  bgmCompare_interface.cpp \
+  bgm_interface.cpp \
+  data_simulation.cpp \
+  math/custom_exp.cpp \
+  math/e_exp.cpp \
+  mcmc/mcmc_hmc.cpp \
+  mcmc/mcmc_leapfrog.cpp \
+  mcmc/mcmc_nuts.cpp \
+  mcmc/mcmc_rwm.cpp \
+  mcmc/mcmc_utils.cpp \
+  priors/sbm_edge_prior.cpp \
+  sbm_edge_prior_interface.cpp \
+  utils/progress_manager.cpp \
+
diff --git a/src/common_helpers.h b/src/utils/common_helpers.h
similarity index 100%
rename from src/common_helpers.h
rename to src/utils/common_helpers.h
diff --git a/src/print_mutex.h b/src/utils/print_mutex.h
similarity index 100%
rename from src/print_mutex.h
rename to src/utils/print_mutex.h
diff --git a/src/progress_manager.cpp b/src/utils/progress_manager.cpp
similarity index 99%
rename from src/progress_manager.cpp
rename to src/utils/progress_manager.cpp
index 5ad13fbd..83610a69 100644
--- a/src/progress_manager.cpp
+++ b/src/utils/progress_manager.cpp
@@ -1,4 +1,4 @@
-#include "progress_manager.h"
+#include "utils/progress_manager.h"
 
 ProgressManager::ProgressManager(int nChains_, int nIter_, int nWarmup_, int printEvery_, int progress_type_, bool useUnicode_)
   : nChains(nChains_), nIter(nIter_ + nWarmup_), nWarmup(nWarmup_), printEvery(printEvery_),
@@ -398,7 +398,7 @@ void ProgressManager::maybePadToLength(std::string& content) const {
 // };
 
 
-// // [[Rcpp::export]]
+// // [[Rcpp::export]] // if uncommented, must move .cpp file to src/ for Rcpp to compile
 // void runMCMC_parallel(int nChains = 4, int nIter = 100, int nWarmup = 100, int progress_type = 2, bool useUnicode = false,
 //   int delay = 20) {
 
diff --git a/src/progress_manager.h b/src/utils/progress_manager.h
similarity index 100%
rename from src/progress_manager.h
rename to src/utils/progress_manager.h
diff --git a/src/variable_helpers.h b/src/utils/variable_helpers.h
similarity index 99%
rename from src/variable_helpers.h
rename to src/utils/variable_helpers.h
index 79becda6..7940d61c 100644
--- a/src/variable_helpers.h
+++ b/src/utils/variable_helpers.h
@@ -1,5 +1,5 @@
 #include <RcppArmadillo.h>
-#include "explog_switch.h"
+#include "math/explog_switch.h"