diff --git a/src/bbb_exchange/asl_multi_te.py b/src/bbb_exchange/asl_multi_te.py index 5101f64..de0a1c9 100644 --- a/src/bbb_exchange/asl_multi_te.py +++ b/src/bbb_exchange/asl_multi_te.py @@ -4,8 +4,9 @@ from fitting_multi_te import ls_fit_volume_multite, prepare_multite_data from fitting_multi_te import bayesian_fit_volume_multite, create_multite_bayesian_config import json -with open("config.json", "r") as file: - config = json.load(file) +from bbb_exchange.core.config_manager import ConfigManager + +config = ConfigManager("config.json").load() """ Multi-TE ASL Processing Script with Bayesian Fitting diff --git a/src/bbb_exchange/asl_single_te.py b/src/bbb_exchange/asl_single_te.py index 5c73c65..df67242 100644 --- a/src/bbb_exchange/asl_single_te.py +++ b/src/bbb_exchange/asl_single_te.py @@ -22,9 +22,9 @@ - nifti images of the fitted CBF, ATT, ABV and ATT_A data for the extended model - print values in csv file (optional) """ +from bbb_exchange.core.config_manager import ConfigManager -with open("config.json", "r") as file: - config = json.load(file) +config = ConfigManager("config.json").load() def create_parameter_config_from_config(): """ Create parameter configuration from config.py diff --git a/src/bbb_exchange/core/config_manager.py b/src/bbb_exchange/core/config_manager.py new file mode 100644 index 0000000..45b5da3 --- /dev/null +++ b/src/bbb_exchange/core/config_manager.py @@ -0,0 +1,48 @@ +import json +import os +import logging + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + + +class ConfigManager: + """ + Lightweight configuration manager for BBB-ASL. + Centralizes JSON loading and prevents KeyError crashes. + """ + + def __init__(self, config_path="config.json"): + self.config_path = config_path + self.config = {} + + def load(self): + if not os.path.exists(self.config_path): + logger.warning(f"{self.config_path} not found. Using empty configuration.") + return {} + + try: + with open(self.config_path, "r") as f: + self.config = json.load(f) + logger.info(f"Loaded configuration from {self.config_path}") + except json.JSONDecodeError as e: + logger.error(f"Invalid JSON format: {e}") + self.config = {} + except Exception as e: + logger.error(f"Error loading configuration: {e}") + self.config = {} + + return self.config + + def get(self, *keys, default=None): + """ + Safe nested access: + config.get("physiological", "T1", default=1.6) + """ + value = self.config + for key in keys: + if isinstance(value, dict) and key in value: + value = value[key] + else: + return default + return value \ No newline at end of file diff --git a/src/bbb_exchange/fitting_multi_te.py b/src/bbb_exchange/fitting_multi_te.py index dfb5720..f7d537f 100644 --- a/src/bbb_exchange/fitting_multi_te.py +++ b/src/bbb_exchange/fitting_multi_te.py @@ -4,6 +4,21 @@ from model_multi_te import deltaM_multite_model import stan +import os + +def load_stan_model(filename): + base_dir = os.path.dirname(os.path.abspath(__file__)) + stan_path = os.path.join(base_dir, "models", "stan", filename) + + if not os.path.exists(stan_path): + raise FileNotFoundError(f"Stan file not found: {stan_path}") + + with open(stan_path, "r") as f: + return f.read() + + +STAN_MODEL_CODE = load_stan_model("multite_full.stan") + import json with open("config.json", "r") as file: config = json.load(file) @@ -40,245 +55,6 @@ def bayesian_fit_voxel_multite(tis, tes, ntes, signal, M0a, taus, cbf_prior_std = 20.0 - - STAN_MODEL_CODE = """ -// multi_te_full.stan -data { - // Measurements - int n_measurements; // total number of TE points over all TI - int n_ti; // number different TI - array[n_measurements] real signal; // measured DeltaM (normalised) - - // Timing - array[n_ti] real tis; // TI[j] - array[n_ti] int ntes; // number TE for each TI - array[n_measurements] real tes; - array[n_ti] real tau_per_ti; // labelling time - - real t1; - real t1b; - real t2; - real t2b; - real texch; - real itt; - real lambd; - real alpha; - real M0a; - - // Priors - real att_prior_mean; - real att_prior_std; - real cbf_prior_mean; - real cbf_prior_std; -} - -parameters { - real att; // s - real cbf; // ml/min/100g - real sigma; // noise -} - -transformed parameters { - // conversion: cbf [ml/min/100g] -> f [ml/s/g] - real f = (cbf / 100.0) * 60.0 / 6000.0; - vector[n_measurements] mu; - - { - int te_index = 1; // Stan-index starts at 1 (not 0) - vector[n_measurements] S_bl1_final = rep_vector(0.0, n_measurements); - vector[n_measurements] S_bl2_final = rep_vector(0.0, n_measurements); - vector[n_measurements] S_ex_final = rep_vector(0.0, n_measurements); - - for (j in 1:n_ti) { - real tau = tau_per_ti[j]; - real ti = tis[j]; - - // === Case 1: 0 < ti < att === - if ((0 < ti) && (ti < att)) { - for (k in 1:ntes[j]) { - S_bl1_final[te_index] = 0; - S_bl2_final[te_index] = 0; - S_ex_final[te_index] = 0; - te_index += 1; - } - } - - // === Case 2: att <= ti < (att + itt) === - else if ((att <= ti) && (ti < att + itt)) { - for (k in 1:ntes[j]) { - real te = tes[te_index]; - if ((0 <= te) && (te < (att + itt - ti))) { - S_bl1_final[te_index] = - (2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) - * (exp(ti / t1b) - exp(att / t1b)) * exp(-te / t2b)); - } - else if (((att + itt - ti) <= te) && (te < itt)) { - real base_term = 2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) - * (exp(ti / t1b) - exp(att / t1b)); - real transition_factor = (te - (att + itt - ti)) / (ti - att); - S_bl1_final[te_index] = (base_term - transition_factor * base_term) * exp(-te / t2b); - S_bl2_final[te_index] = (transition_factor * base_term) * exp(-te / t2b) * exp(-te / texch); - S_ex_final[te_index] = (transition_factor * base_term) * (1 - exp(-te / texch)) * exp(-te / t2); - } - else { - S_bl2_final[te_index] = - (2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) - * (exp(ti / t1b) - exp(att / t1b)) * exp(-te / t2b) * exp(-te / texch)); - S_ex_final[te_index] = - (2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) - * (exp(ti / t1b) - exp(att / t1b)) * (1 - exp(-te / texch)) * exp(-te / t2)); - } - te_index += 1; - } - } - - // === Case 3: (att+itt) <= ti < (att + tau) === - else if (((att + itt) <= ti) && (ti < (att + tau))) { - for (k in 1:ntes[j]) { - real te = tes[te_index]; - real term1 = 2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) - * (exp(ti / t1b) - exp(att / t1b)); - real term2 = 2 * f * t1b * exp(-(att + itt) / t1b) * exp(-ti / t1b) - * (exp(ti / t1b) - exp((att + itt) / t1b)); - real base_diff = term1 - term2; - if ((0 <= te) && (te < itt)) { - real transition_factor = te / itt; - S_bl1_final[te_index] = (base_diff - transition_factor * base_diff) * exp(-te / t2b); - S_bl2_final[te_index] = - ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + transition_factor * base_diff) * exp(-te / t2b) * exp(-te / texch)); - S_ex_final[te_index] = - (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) - * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) - - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) - * exp(-((1 / t1) + (1 / texch)) * ti) - * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) - * exp(-te / t2)) - + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); - } else { - S_bl2_final[te_index] = - ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + base_diff) * exp(-te / t2b) * exp(-te / texch)); - S_ex_final[te_index] = - (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) - * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) - - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) - * exp(-((1 / t1) + (1 / texch)) * ti) - * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) - * exp(-te / t2)) - + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); - } - te_index += 1; - } - } - - // === Case 4: (att + tau) <= ti < (att + itt + tau) === - else if (((att + tau) <= ti) && (ti < (att + itt + tau))) { - for (k in 1:ntes[j]) { - real te = tes[te_index]; - real term1 = 2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) - * (exp((att + tau) / t1b) - exp(att / t1b)); - real term2 = 2 * f * t1b * exp(-(att + itt) / t1b) * exp(-ti / t1b) - * (exp(ti / t1b) - exp((att + itt) / t1b)); - real base_diff = term1 - term2; - if ((0 <= te) && (te < (itt - (ti - (att + tau))))) { - real transition_factor = te / (itt - (ti - (att + tau))); - S_bl1_final[te_index] = (base_diff - transition_factor * base_diff) * exp(-te / t2b); - S_bl2_final[te_index] = - ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + transition_factor * base_diff) * exp(-te / t2b) * exp(-te / texch)); - S_ex_final[te_index] = - (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) - * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) - - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) - * exp(-((1 / t1) + (1 / texch)) * ti) - * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) - * exp(-te / t2)) - + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + transition_factor * base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); - } else { - S_bl2_final[te_index] = - ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + base_diff) * exp(-te / t2b) * exp(-te / texch)); - S_ex_final[te_index] = - (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) - * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) - - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) - * exp(-((1 / t1) + (1 / texch)) * ti) - * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) - * exp(-te / t2)) - + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - + base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); - } - te_index += 1; - } - } - - // === Case 5: ti >= (att + itt + tau) === - else { - for (k in 1:ntes[j]) { - real te = tes[te_index]; - S_bl2_final[te_index] = - (2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * (att + itt + tau)) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) - * exp(-te / t2b) * exp(-te / texch)); - S_ex_final[te_index] = - (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) - * (exp((1 / t1) * (att + itt + tau)) - exp((1 / t1) * (att + itt))) - - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) - * exp(-((1 / t1) + (1 / texch)) * ti) - * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) - * exp(-te / t2)) - + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) - * exp(-((1 / t1b) + (1 / texch)) * ti) - * (exp(((1 / t1b) + (1 / texch)) * (att + itt + tau)) - exp(((1 / t1b) + (1 / texch)) * (att + itt)))) - * (1 - exp(-te / texch)) * exp(-te / t2)); - te_index += 1; - } - } - } - - mu = (S_bl1_final + S_bl2_final + S_ex_final) * (M0a * alpha / lambd); - } -} - -model { - // Priors - att ~ normal(att_prior_mean, att_prior_std); - cbf ~ normal(cbf_prior_mean, cbf_prior_std); - sigma ~ exponential(1.0); - - // Likelihood - signal ~ normal(mu, sigma); -} - -generated quantities { - vector[n_measurements] mu_ppc; - for (i in 1:n_measurements) { - mu_ppc[i] = normal_rng(mu[i], sigma); - } -} - """ - # Prepare data for Stan data = { "n_measurements": len(signal), diff --git a/src/bbb_exchange/fitting_single_te.py b/src/bbb_exchange/fitting_single_te.py index 2a6389e..78bfc6a 100644 --- a/src/bbb_exchange/fitting_single_te.py +++ b/src/bbb_exchange/fitting_single_te.py @@ -5,6 +5,22 @@ from scipy.optimize import curve_fit from DeltaM_model import DeltaM_model_ext, dm_tiss + +def load_stan_model(filename): + base_dir = os.path.dirname(os.path.abspath(__file__)) + stan_path = os.path.join(base_dir, "models", "stan", filename) + + if not os.path.exists(stan_path): + raise FileNotFoundError(f"Stan file not found: {stan_path}") + + with open(stan_path, "r") as f: + return f.read() + + +STAN_MODEL_CODE = load_stan_model("single_compartment.stan") +STAN_MODEL_CODE_EXT = load_stan_model("extended_compartment.stan") + + import json """ @@ -218,137 +234,6 @@ def ls_fit_volume_ext(pwi_data, t, m0_data, tau, lambd=0.9, T1=1.6, T1a=1.65, a= # === Fitting functions for Bayesian fitting -STAN_MODEL_CODE = """ -functions { - vector deltaM_model(vector t, real att, real f, real M0a, real tau, - real T1a, real T1, real lambd, real a) { - int n = num_elements(t); - vector[n] deltaM = rep_vector(0.0, n); - - real T1app = 1.0 / (1.0 / T1 + f / lambd); - real R = 1.0 / T1app - 1.0 / T1a; - - for (i in 1:n) { - if (t[i] >= att && t[i] <= att + tau) { - // Case 2: att <= t <= att + tau - real time = t[i]; - real term = exp(R * time) - exp(R * att); - deltaM[i] = (2 * a * M0a * f * exp(-time / T1app) / R) * term; - } else if (t[i] > att + tau) { - // Case 3: t > att + tau - real time = t[i]; - real term = exp(R * (att + tau)) - exp(R * att); - deltaM[i] = (2 * a * M0a * f * exp(-time / T1app) / R) * term; - } - // else: t < att, deltaM[i] remains 0 - } - return deltaM; - } -} - -data { - int n; - vector[n] t; - vector[n] signal; - real M0a; - real tau; - - // Fixed parameter values - real T1a_fixed; - real T1_fixed; - real lambd_fixed; - real a_fixed; - - // Fitting flags - int fit_T1a; - int fit_T1; - int fit_lambd; - - // Priors - real T1a_prior_mean; - real T1a_prior_std; - real T1_prior_mean; - real T1_prior_std; - real lambd_prior_mean; - real lambd_prior_std; - - // ATT and CBF priors from LS fitting - int use_att_prior_from_ls; - int use_cbf_prior_from_ls; - real att_prior_from_ls; - real att_prior_std; - real cbf_prior_from_ls; - real cbf_prior_std; - - // Bounds - real T1a_lower; - real T1a_upper; - real T1_lower; - real T1_upper; - real lambd_lower; - real lambd_upper; -} - -parameters { - real att; - real f; // CBF as f parameter - real sigma; - - // Conditional parameters - real T1a_param; - real T1_param; - real lambd_param; -} - -transformed parameters { - real T1a_use; - real T1_use; - real lambd_use; - real cbf; // CBF in ml/min/100g - - // Use fitted or fixed values - T1a_use = fit_T1a ? T1a_param : T1a_fixed; - T1_use = fit_T1 ? T1_param : T1_fixed; - lambd_use = fit_lambd ? lambd_param : lambd_fixed; - - // Convert f to CBF - cbf = f * 6000.0; -} - -model { - vector[n] mu; - - if (use_att_prior_from_ls == 1) { - att ~ normal(att_prior_from_ls, att_prior_std); - } else { - att ~ normal(1.2, 0.5); - } - - if (use_cbf_prior_from_ls == 1) { - f ~ normal(cbf_prior_from_ls / 6000.0, cbf_prior_std / 6000.0); - } else { - f ~ normal(0.01, 0.0025); // Default: 60 ml/min/100g ± 15 - } - - sigma ~ exponential(1); - - // Conditional priors for flexible parameters - if (fit_T1a) { - T1a_param ~ normal(T1a_prior_mean, T1a_prior_std); - } - if (fit_T1) { - T1_param ~ normal(T1_prior_mean, T1_prior_std); - } - if (fit_lambd) { - lambd_param ~ normal(lambd_prior_mean, lambd_prior_std); - } - - // Likelihood - mu = deltaM_model(t, att, f, M0a, tau, T1a_use, T1_use, lambd_use, a_fixed); - signal ~ normal(mu, sigma); -} -""" - def bayesian_fit_voxel(t, signal, M0a, tau, T1a_val, T1_val, lambd_val, a_val, @@ -582,166 +467,6 @@ def bayesian_fit_volume(pwi_data, t, m0_data, tau, param_config, att_ls_map=None - -STAN_MODEL_CODE_EXT = """ -functions { - vector deltaM_model_ext(vector t, real att, real f, real M0a, real tau, - real abv, real att_a, real T1a, real T1, real lambd, real a) { - int n = num_elements(t); - vector[n] deltaM = rep_vector(0.0, n); - - // Tissue compartment (Eq. 1) - real T1app = 1.0 / (1.0 / T1 + f / lambd); - real R = 1.0 / T1app - 1.0 / T1a; - - for (i in 1:n) { - // Tissue - if (t[i] >= att && t[i] <= att + tau) { - real term = exp(R * t[i]) - exp(R * att); - deltaM[i] += (2 * a * M0a * f * exp(-t[i] / T1app) / R) * term; - } else if (t[i] > att + tau) { - real term = exp(R * (att + tau)) - exp(R * att); - deltaM[i] += (2 * a * M0a * f * exp(-t[i] / T1app) / R) * term; - } - - // Arterial compartment (Eq. 2) - if (t[i] >= att_a && t[i] <= att_a + tau) { - deltaM[i] += 2 * a * M0a * abv * exp(-t[i] / T1a); - } - } - return deltaM; - } -} - -data { - int n; - vector[n] t; - vector[n] signal; - real M0a; - real tau; - - // Parameter values - real T1a_fixed; - real T1_fixed; - real lambd_fixed; - real a_fixed; - real abv_fixed; - real att_a_fixed; - - // Fitting flags - int fit_T1a; - int fit_T1; - int fit_lambd; - int fit_abv; - int fit_att_a; - - // Priors - real T1a_prior_mean; - real T1a_prior_std; - real T1_prior_mean; - real T1_prior_std; - real lambd_prior_mean; - real lambd_prior_std; - real abv_prior_mean; - real abv_prior_std; - real att_a_prior_mean; - real att_a_prior_std; - - // ATT and CBF priors from LS fitting - int use_att_prior_from_ls; - int use_cbf_prior_from_ls; - real att_prior_from_ls; - real att_prior_std; - real cbf_prior_from_ls; - real cbf_prior_std; - - // Bounds - real T1a_lower; - real T1a_upper; - real T1_lower; - real T1_upper; - real lambd_lower; - real lambd_upper; - real abv_lower; - real abv_upper; - real att_a_lower; - real att_a_upper; -} - -parameters { - real att; - real f; // CBF as f parameter - real sigma; - - // Conditional parameters - real T1a_param; - real T1_param; - real lambd_param; - real abv_param; - real att_a_param; -} - -transformed parameters { - real T1a_use; - real T1_use; - real lambd_use; - real abv_use; - real att_a_use; - real cbf; // CBF in ml/min/100g - - // Use fitted or constant fixed values - T1a_use = fit_T1a ? T1a_param : T1a_fixed; - T1_use = fit_T1 ? T1_param : T1_fixed; - lambd_use = fit_lambd ? lambd_param : lambd_fixed; - abv_use = fit_abv ? abv_param : abv_fixed; - att_a_use = fit_att_a ? att_a_param : att_a_fixed; - - // Convert f to CBF - cbf = f * 6000.0; -} - -model { - vector[n] mu; - - // Priors for ATT and CBF - if (use_att_prior_from_ls == 1) { - att ~ normal(att_prior_from_ls, att_prior_std); - } else { - att ~ normal(1.2, 0.5); - } - - if (use_cbf_prior_from_ls == 1) { - f ~ normal(cbf_prior_from_ls / 6000.0, cbf_prior_std / 6000.0); - } else { - f ~ normal(0.01, 0.0025); // Default: 60 ml/min/100g - } - - sigma ~ exponential(1); - - // Conditional priors for flexible parameters - if (fit_T1a) { - T1a_param ~ normal(T1a_prior_mean, T1a_prior_std); - } - if (fit_T1) { - T1_param ~ normal(T1_prior_mean, T1_prior_std); - } - if (fit_lambd) { - lambd_param ~ normal(lambd_prior_mean, lambd_prior_std); - } - if (fit_abv) { - abv_param ~ normal(abv_prior_mean, abv_prior_std); - } - if (fit_att_a) { - att_a_param ~ normal(att_a_prior_mean, att_a_prior_std); - } - - // Likelihood - mu = deltaM_model_ext(t, att, f, M0a, tau, abv_use, att_a_use, - T1a_use, T1_use, lambd_use, a_fixed); - signal ~ normal(mu, sigma); -} -""" - def bayesian_fit_voxel_ext(t, signal, M0a, tau, T1a_val, T1_val, lambd_val, a_val, param_config, diff --git a/src/bbb_exchange/models/stan/extended_compartment.stan b/src/bbb_exchange/models/stan/extended_compartment.stan new file mode 100644 index 0000000..e6ea242 --- /dev/null +++ b/src/bbb_exchange/models/stan/extended_compartment.stan @@ -0,0 +1,156 @@ +functions { + vector deltaM_model_ext(vector t, real att, real f, real M0a, real tau, + real abv, real att_a, real T1a, real T1, real lambd, real a) { + int n = num_elements(t); + vector[n] deltaM = rep_vector(0.0, n); + + // Tissue compartment (Eq. 1) + real T1app = 1.0 / (1.0 / T1 + f / lambd); + real R = 1.0 / T1app - 1.0 / T1a; + + for (i in 1:n) { + // Tissue + if (t[i] >= att && t[i] <= att + tau) { + real term = exp(R * t[i]) - exp(R * att); + deltaM[i] += (2 * a * M0a * f * exp(-t[i] / T1app) / R) * term; + } else if (t[i] > att + tau) { + real term = exp(R * (att + tau)) - exp(R * att); + deltaM[i] += (2 * a * M0a * f * exp(-t[i] / T1app) / R) * term; + } + + // Arterial compartment (Eq. 2) + if (t[i] >= att_a && t[i] <= att_a + tau) { + deltaM[i] += 2 * a * M0a * abv * exp(-t[i] / T1a); + } + } + return deltaM; + } +} + +data { + int n; + vector[n] t; + vector[n] signal; + real M0a; + real tau; + + // Parameter values + real T1a_fixed; + real T1_fixed; + real lambd_fixed; + real a_fixed; + real abv_fixed; + real att_a_fixed; + + // Fitting flags + int fit_T1a; + int fit_T1; + int fit_lambd; + int fit_abv; + int fit_att_a; + + // Priors + real T1a_prior_mean; + real T1a_prior_std; + real T1_prior_mean; + real T1_prior_std; + real lambd_prior_mean; + real lambd_prior_std; + real abv_prior_mean; + real abv_prior_std; + real att_a_prior_mean; + real att_a_prior_std; + + // ATT and CBF priors from LS fitting + int use_att_prior_from_ls; + int use_cbf_prior_from_ls; + real att_prior_from_ls; + real att_prior_std; + real cbf_prior_from_ls; + real cbf_prior_std; + + // Bounds + real T1a_lower; + real T1a_upper; + real T1_lower; + real T1_upper; + real lambd_lower; + real lambd_upper; + real abv_lower; + real abv_upper; + real att_a_lower; + real att_a_upper; +} + +parameters { + real att; + real f; // CBF as f parameter + real sigma; + + // Conditional parameters + real T1a_param; + real T1_param; + real lambd_param; + real abv_param; + real att_a_param; +} + +transformed parameters { + real T1a_use; + real T1_use; + real lambd_use; + real abv_use; + real att_a_use; + real cbf; // CBF in ml/min/100g + + // Use fitted or constant fixed values + T1a_use = fit_T1a ? T1a_param : T1a_fixed; + T1_use = fit_T1 ? T1_param : T1_fixed; + lambd_use = fit_lambd ? lambd_param : lambd_fixed; + abv_use = fit_abv ? abv_param : abv_fixed; + att_a_use = fit_att_a ? att_a_param : att_a_fixed; + + // Convert f to CBF + cbf = f * 6000.0; +} + +model { + vector[n] mu; + + // Priors for ATT and CBF + if (use_att_prior_from_ls == 1) { + att ~ normal(att_prior_from_ls, att_prior_std); + } else { + att ~ normal(1.2, 0.5); + } + + if (use_cbf_prior_from_ls == 1) { + f ~ normal(cbf_prior_from_ls / 6000.0, cbf_prior_std / 6000.0); + } else { + f ~ normal(0.01, 0.0025); // Default: 60 ml/min/100g + } + + sigma ~ exponential(1); + + // Conditional priors for flexible parameters + if (fit_T1a) { + T1a_param ~ normal(T1a_prior_mean, T1a_prior_std); + } + if (fit_T1) { + T1_param ~ normal(T1_prior_mean, T1_prior_std); + } + if (fit_lambd) { + lambd_param ~ normal(lambd_prior_mean, lambd_prior_std); + } + if (fit_abv) { + abv_param ~ normal(abv_prior_mean, abv_prior_std); + } + if (fit_att_a) { + att_a_param ~ normal(att_a_prior_mean, att_a_prior_std); + } + + // Likelihood + mu = deltaM_model_ext(t, att, f, M0a, tau, abv_use, att_a_use, + T1a_use, T1_use, lambd_use, a_fixed); + signal ~ normal(mu, sigma); +} \ No newline at end of file diff --git a/src/bbb_exchange/models/stan/multite_full.stan b/src/bbb_exchange/models/stan/multite_full.stan new file mode 100644 index 0000000..f66d55b --- /dev/null +++ b/src/bbb_exchange/models/stan/multite_full.stan @@ -0,0 +1,235 @@ +// multi_te_full.stan +data { + // Measurements + int n_measurements; // total number of TE points over all TI + int n_ti; // number different TI + array[n_measurements] real signal; // measured DeltaM (normalised) + + // Timing + array[n_ti] real tis; // TI[j] + array[n_ti] int ntes; // number TE for each TI + array[n_measurements] real tes; + array[n_ti] real tau_per_ti; // labelling time + + real t1; + real t1b; + real t2; + real t2b; + real texch; + real itt; + real lambd; + real alpha; + real M0a; + + // Priors + real att_prior_mean; + real att_prior_std; + real cbf_prior_mean; + real cbf_prior_std; +} + +parameters { + real att; // s + real cbf; // ml/min/100g + real sigma; // noise +} + +transformed parameters { + // conversion: cbf [ml/min/100g] -> f [ml/s/g] + real f = (cbf / 100.0) * 60.0 / 6000.0; + vector[n_measurements] mu; + + { + int te_index = 1; // Stan-index starts at 1 (not 0) + vector[n_measurements] S_bl1_final = rep_vector(0.0, n_measurements); + vector[n_measurements] S_bl2_final = rep_vector(0.0, n_measurements); + vector[n_measurements] S_ex_final = rep_vector(0.0, n_measurements); + + for (j in 1:n_ti) { + real tau = tau_per_ti[j]; + real ti = tis[j]; + + // === Case 1: 0 < ti < att === + if ((0 < ti) && (ti < att)) { + for (k in 1:ntes[j]) { + S_bl1_final[te_index] = 0; + S_bl2_final[te_index] = 0; + S_ex_final[te_index] = 0; + te_index += 1; + } + } + + // === Case 2: att <= ti < (att + itt) === + else if ((att <= ti) && (ti < att + itt)) { + for (k in 1:ntes[j]) { + real te = tes[te_index]; + if ((0 <= te) && (te < (att + itt - ti))) { + S_bl1_final[te_index] = + (2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) + * (exp(ti / t1b) - exp(att / t1b)) * exp(-te / t2b)); + } + else if (((att + itt - ti) <= te) && (te < itt)) { + real base_term = 2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) + * (exp(ti / t1b) - exp(att / t1b)); + real transition_factor = (te - (att + itt - ti)) / (ti - att); + S_bl1_final[te_index] = (base_term - transition_factor * base_term) * exp(-te / t2b); + S_bl2_final[te_index] = (transition_factor * base_term) * exp(-te / t2b) * exp(-te / texch); + S_ex_final[te_index] = (transition_factor * base_term) * (1 - exp(-te / texch)) * exp(-te / t2); + } + else { + S_bl2_final[te_index] = + (2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) + * (exp(ti / t1b) - exp(att / t1b)) * exp(-te / t2b) * exp(-te / texch)); + S_ex_final[te_index] = + (2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) + * (exp(ti / t1b) - exp(att / t1b)) * (1 - exp(-te / texch)) * exp(-te / t2)); + } + te_index += 1; + } + } + + // === Case 3: (att+itt) <= ti < (att + tau) === + else if (((att + itt) <= ti) && (ti < (att + tau))) { + for (k in 1:ntes[j]) { + real te = tes[te_index]; + real term1 = 2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) + * (exp(ti / t1b) - exp(att / t1b)); + real term2 = 2 * f * t1b * exp(-(att + itt) / t1b) * exp(-ti / t1b) + * (exp(ti / t1b) - exp((att + itt) / t1b)); + real base_diff = term1 - term2; + if ((0 <= te) && (te < itt)) { + real transition_factor = te / itt; + S_bl1_final[te_index] = (base_diff - transition_factor * base_diff) * exp(-te / t2b); + S_bl2_final[te_index] = + ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + transition_factor * base_diff) * exp(-te / t2b) * exp(-te / texch)); + S_ex_final[te_index] = + (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) + * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) + - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) + * exp(-((1 / t1) + (1 / texch)) * ti) + * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) + * exp(-te / t2)) + + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); + } else { + S_bl2_final[te_index] = + ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + base_diff) * exp(-te / t2b) * exp(-te / texch)); + S_ex_final[te_index] = + (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) + * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) + - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) + * exp(-((1 / t1) + (1 / texch)) * ti) + * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) + * exp(-te / t2)) + + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); + } + te_index += 1; + } + } + + // === Case 4: (att + tau) <= ti < (att + itt + tau) === + else if (((att + tau) <= ti) && (ti < (att + itt + tau))) { + for (k in 1:ntes[j]) { + real te = tes[te_index]; + real term1 = 2 * f * t1b * exp(-att / t1b) * exp(-ti / t1b) + * (exp((att + tau) / t1b) - exp(att / t1b)); + real term2 = 2 * f * t1b * exp(-(att + itt) / t1b) * exp(-ti / t1b) + * (exp(ti / t1b) - exp((att + itt) / t1b)); + real base_diff = term1 - term2; + if ((0 <= te) && (te < (itt - (ti - (att + tau))))) { + real transition_factor = te / (itt - (ti - (att + tau))); + S_bl1_final[te_index] = (base_diff - transition_factor * base_diff) * exp(-te / t2b); + S_bl2_final[te_index] = + ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + transition_factor * base_diff) * exp(-te / t2b) * exp(-te / texch)); + S_ex_final[te_index] = + (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) + * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) + - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) + * exp(-((1 / t1) + (1 / texch)) * ti) + * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) + * exp(-te / t2)) + + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + transition_factor * base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); + } else { + S_bl2_final[te_index] = + ((2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + base_diff) * exp(-te / t2b) * exp(-te / texch)); + S_ex_final[te_index] = + (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) + * (exp((1 / t1) * ti) - exp((1 / t1) * (att + itt))) + - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) + * exp(-((1 / t1) + (1 / texch)) * ti) + * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) + * exp(-te / t2)) + + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * ti) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + + base_diff) * (1 - exp(-te / texch)) * exp(-te / t2)); + } + te_index += 1; + } + } + + // === Case 5: ti >= (att + itt + tau) === + else { + for (k in 1:ntes[j]) { + real te = tes[te_index]; + S_bl2_final[te_index] = + (2 * f * exp(-(1 / t1b) * (att + itt)) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * (att + itt + tau)) - exp(((1 / t1b) + (1 / texch)) * (att + itt))) + * exp(-te / t2b) * exp(-te / texch)); + S_ex_final[te_index] = + (((2 * f * exp(-(1 / t1b) * (att + itt))) / (1 / t1) * exp(-(1 / t1) * ti) + * (exp((1 / t1) * (att + itt + tau)) - exp((1 / t1) * (att + itt))) + - (2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / texch) + (1 / t1)) + * exp(-((1 / t1) + (1 / texch)) * ti) + * (exp(((1 / texch) + (1 / t1)) * ti) - exp(((1 / texch) + (1 / t1)) * (att + itt)))) + * exp(-te / t2)) + + (((2 * f * exp(-(1 / t1b) * (att + itt))) / ((1 / t1b) + (1 / texch)) + * exp(-((1 / t1b) + (1 / texch)) * ti) + * (exp(((1 / t1b) + (1 / texch)) * (att + itt + tau)) - exp(((1 / t1b) + (1 / texch)) * (att + itt)))) + * (1 - exp(-te / texch)) * exp(-te / t2)); + te_index += 1; + } + } + } + + mu = (S_bl1_final + S_bl2_final + S_ex_final) * (M0a * alpha / lambd); + } +} + +model { + // Priors + att ~ normal(att_prior_mean, att_prior_std); + cbf ~ normal(cbf_prior_mean, cbf_prior_std); + sigma ~ exponential(1.0); + + // Likelihood + signal ~ normal(mu, sigma); +} + +generated quantities { + vector[n_measurements] mu_ppc; + for (i in 1:n_measurements) { + mu_ppc[i] = normal_rng(mu[i], sigma); + } +} \ No newline at end of file diff --git a/src/bbb_exchange/models/stan/single_compartment.stan b/src/bbb_exchange/models/stan/single_compartment.stan new file mode 100644 index 0000000..1c17832 --- /dev/null +++ b/src/bbb_exchange/models/stan/single_compartment.stan @@ -0,0 +1,128 @@ +functions { + vector deltaM_model(vector t, real att, real f, real M0a, real tau, + real T1a, real T1, real lambd, real a) { + int n = num_elements(t); + vector[n] deltaM = rep_vector(0.0, n); + + real T1app = 1.0 / (1.0 / T1 + f / lambd); + real R = 1.0 / T1app - 1.0 / T1a; + + for (i in 1:n) { + if (t[i] >= att && t[i] <= att + tau) { + // Case 2: att <= t <= att + tau + real time = t[i]; + real term = exp(R * time) - exp(R * att); + deltaM[i] = (2 * a * M0a * f * exp(-time / T1app) / R) * term; + } else if (t[i] > att + tau) { + // Case 3: t > att + tau + real time = t[i]; + real term = exp(R * (att + tau)) - exp(R * att); + deltaM[i] = (2 * a * M0a * f * exp(-time / T1app) / R) * term; + } + // else: t < att, deltaM[i] remains 0 + } + return deltaM; + } +} + +data { + int n; + vector[n] t; + vector[n] signal; + real M0a; + real tau; + + // Fixed parameter values + real T1a_fixed; + real T1_fixed; + real lambd_fixed; + real a_fixed; + + // Fitting flags + int fit_T1a; + int fit_T1; + int fit_lambd; + + // Priors + real T1a_prior_mean; + real T1a_prior_std; + real T1_prior_mean; + real T1_prior_std; + real lambd_prior_mean; + real lambd_prior_std; + + // ATT and CBF priors from LS fitting + int use_att_prior_from_ls; + int use_cbf_prior_from_ls; + real att_prior_from_ls; + real att_prior_std; + real cbf_prior_from_ls; + real cbf_prior_std; + + // Bounds + real T1a_lower; + real T1a_upper; + real T1_lower; + real T1_upper; + real lambd_lower; + real lambd_upper; +} + +parameters { + real att; + real f; // CBF as f parameter + real sigma; + + // Conditional parameters + real T1a_param; + real T1_param; + real lambd_param; +} + +transformed parameters { + real T1a_use; + real T1_use; + real lambd_use; + real cbf; // CBF in ml/min/100g + + // Use fitted or fixed values + T1a_use = fit_T1a ? T1a_param : T1a_fixed; + T1_use = fit_T1 ? T1_param : T1_fixed; + lambd_use = fit_lambd ? lambd_param : lambd_fixed; + + // Convert f to CBF + cbf = f * 6000.0; +} + +model { + vector[n] mu; + + if (use_att_prior_from_ls == 1) { + att ~ normal(att_prior_from_ls, att_prior_std); + } else { + att ~ normal(1.2, 0.5); + } + + if (use_cbf_prior_from_ls == 1) { + f ~ normal(cbf_prior_from_ls / 6000.0, cbf_prior_std / 6000.0); + } else { + f ~ normal(0.01, 0.0025); // Default: 60 ml/min/100g ± 15 + } + + sigma ~ exponential(1); + + // Conditional priors for flexible parameters + if (fit_T1a) { + T1a_param ~ normal(T1a_prior_mean, T1a_prior_std); + } + if (fit_T1) { + T1_param ~ normal(T1_prior_mean, T1_prior_std); + } + if (fit_lambd) { + lambd_param ~ normal(lambd_prior_mean, lambd_prior_std); + } + + // Likelihood + mu = deltaM_model(t, att, f, M0a, tau, T1a_use, T1_use, lambd_use, a_fixed); + signal ~ normal(mu, sigma); +} \ No newline at end of file diff --git a/src/bbb_exchange/stan_loader.py b/src/bbb_exchange/stan_loader.py new file mode 100644 index 0000000..fae50d9 --- /dev/null +++ b/src/bbb_exchange/stan_loader.py @@ -0,0 +1,15 @@ +import os + + +def load_stan_model(stan_filename): + """ + Load Stan model code from file. + """ + base_dir = os.path.dirname(os.path.abspath(__file__)) + stan_path = os.path.join(base_dir, "models", "stan", stan_filename) + + if not os.path.exists(stan_path): + raise FileNotFoundError(f"Stan model file not found: {stan_path}") + + with open(stan_path, "r") as f: + return f.read() \ No newline at end of file