OSCR

Axon Diameter Mapping in the Living Human Brain with Ultra-High-Gradient Diffusion MRI at 500 mT/m Gradient Strength.

Code ↔ Paper

5 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 5 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › AxCaliber‐SMT Model Fitting › Forward Model and Assumptions ↔ AxCaliberSMT/gpuAxCaliberSMT.m, lines 81–127 · score 0.83 · AxCaliber, SMT model, axial diffusivity, CSF diffusivity, intrinsic diffusivity, extra cellular
  2. [2] § Methods › AxCaliber‐SMT Model Fitting › Markov Chain Monte Carlo With Ensemble Sampler With Affine Invariance ↔ utils/mcmc.m, lines 385–442 · score 0.64 · affine invariant ensemble, sampler, Goodman, Weare, MCMC, dimensional
  3. [3] § Methods › AxCaliber‐SMT Model Fitting › Markov Chain Monte Carlo With Ensemble Sampler With Affine Invariance ↔ utils/mcmc.m, lines 385–442 · score 0.57 · stretch move, walkers, thinned, Sampler, Affine, Ensemble
  4. [4] § Methods › MRI Acquisition ↔ helpers/diffusion_parameters_exporter.py, the whole file · a weak match · score 0.53 · encoding directions, pulse width, 30 ms, volumes, diffusion
  5. [5] § Methods › Noise Propagation From Synthetic Diffusion Signal ↔ MCRMWI/gpuMCRMWI.m, lines 566–714 · score 0.51 · intra axonal, volume fraction, forward model, component, extracellular, water

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 878 lines · 48 KB · GPL-3.0 · 2 matches

  1. classdef mcmc < handle
  2. % Kwok-Shing Chan @ MGH
  3. % [email hidden]
  4. %
  5. % This is the class of all MCMC related functions
  6. %
  7. % Date created: 13 June 2024
  8. % Date modified: 7 August 2024
  9. % Date modified: 23 August 2024
  10. % Date modified: 5 October 2024
  11. % Date modified: 4 June 2026 (update affine-invariant method with and without global parameters)
  12. %
  13. properties (GetAccess = public, SetAccess = protected)
  15. end
  17. methods
  18. function out = optimisation(this, data, mask, weights, pars0, fitting, FWDfunc, varargin)
  19. % Input
  20. % ----------
  21. % data : N-D measurement data, First 3 dims reserve for spatial info
  22. % mask : M-D signal mask (M=[1,3])
  23. % weights : N-D weights for optimisaiton, same dim as 'data'
  24. % pars0 : Structure variable containing all parameters to be estimated
  25. % fitting : Structure variable containing all fitting algorithm setting
  26. % .modelParams : 1xM cell variable, name of the model parameters, e.g. {'S0','R2star','noise'};
  27. % .lb : 1xM numeric variable, fitting lower bound, same order as field 'modelParams', e.g. [0.5, 0, 0.001];
  28. % .ub : 1xM numeric variable, fitting upper bound, same order as field 'modelParams', e.g. [2, 1, 0.1];
  29. % .algorithm : MCMC algorithm, 'MH'|'GW'
  30. % .iteration : # MCMC iterations
  31. % .thinning : sampling interval between iterations
  32. % .burnin : iterations at the beginning to be discarded, if burnin>1, then the exact number will be used; if 0<burnin<1 then actual burnin = iteration*burnin
  33. % .repetition : # repetition of MCMC proposal
  34. % .xStepSize : step size of model parameter in MCMC proposal, same size and order as 'modelParams' ('MH' only)
  35. % .StepSize : step size for 'GW' in MCMC proposal ('GW' only)
  36. % .Nwalker : # random walkers ('GW' only)
  37. % FWDfunc : function handle of forward model
  38. % varargin : contains additional input requires for FWDfunc
  39. %
  41. fitting = this.check_set_default_basic(fitting);
  43. % Step 0: display basic messages
  44. this.display_basic_algorithm_parameters(fitting);
  46. % mask data to reduce memory load
  47. mask_idx = find(mask>0);
  48. if ~ismatrix(data); data = utils.reshape_ND2GD(data, mask_idx); else; data = data(:,mask_idx); end
  49. if ~ismatrix(weights); weights = utils.reshape_ND2GD(weights, mask_idx); elseif ~isempty(weights); weights = weights(:,mask_idx); end
  50. % data = utils.reshape_ND2GD(data,mask);
  51. % if ~isempty(weights); weights = utils.reshape_ND2GD(weights,mask); else; weights = ones(size(data), 'like', data); end
  52. pars0 = utils.reshape_ND2GD_struct(pars0,mask);
  54. isGlobal = false;
  55. for k = 1:numel(fitting.modelParams)
  56. if k == 1
  57. N = numel(pars0.(fitting.modelParams{k}));
  58. else
  59. if N ~= numel(pars0.(fitting.modelParams{k})) && numel(pars0.(fitting.modelParams{k})) == 1
  60. isGlobal = true;
  61. break;
  62. end
  63. end
  64. end
  66. % MCMC
  67. if strcmpi(fitting.algorithm,'mh')
  68. xPosterior = this.metropolis_hastings(data, pars0, weights, fitting, FWDfunc ,varargin{:});
  69. else
  70. if ~isGlobal
  71. xPosterior = this.goodman_weare(data, pars0, weights, fitting, FWDfunc ,varargin{:});
  72. else
  73. xPosterior = this.goodman_weare_wglobal_constant(data, pars0, weights, fitting, FWDfunc ,varargin{:});
  74. end
  75. end
  77. % finish up
  78. out = this.res2out(xPosterior,fitting,mask);
  80. end
  82. function xPosterior = metropolis_hastings(this,y,x0,weights,fitting,FWDfunc,varargin)
  83. % Input
  84. % ------
  85. % y : measurements, [Nmeas,Nvoxels]
  86. % x0 : structure array, starting points, N fields, each field 1xNvoxel
  87. % weights : weighting for non-linear least square fitting, same dimension as y
  88. % fitting : Structure variable containing all fitting algorithm setting
  89. % .modelParams : 1xM cell variable, name of the model parameters, e.g. {'S0','R2star','noise'};
  90. % .lb : 1xM numeric variable, fitting lower bound, same order as field 'modelParams', e.g. [0.5, 0, 0.001];
  91. % .ub : 1xM numeric variable, fitting upper bound, same order as field 'modelParams', e.g. [2, 1, 0.1];
  92. % .iteration : # MCMC iterations
  93. % .thinning : sampling interval between iterations
  94. % .burnin : iterations at the beginning to be discarded, if burnin>1, then the exact number will be used; if 0<burnin<1 then actual burnin = iteration*burnin
  95. % .repetition : # repetition of MCMC proposal
  96. % .xStepSize : step size of model parameter in MCMC proposal, same size and order as 'modelParams' ('MH' only)
  97. % FWDfunc : function handle for forward signal model
  98. % varargin : other input required for @FWDfunc
  99. %
  101. fitting = this.check_set_default_basic(fitting);
  102. if isempty(weights); weights = ones(size(y), 'like', y); end
  104. % Nm: # measurements; Nv: # voxels
  105. [Nm, Nv] = size(y);
  106. % Nvar: # estimation parameters
  107. Nvar = numel(fitting.modelParams);
  108. Nburnin = this.get_number_burnin(fitting);
  109. % Ns: # samples in posterior distribution
  110. Ns = numel(Nburnin+1:fitting.thinning:fitting.iteration); %floor( (fitting.iteration - floor(fitting.iteration*fitting.burnin)) / fitting.thinning );
  112. % convert data into single datatype for better performance and out themn into GPU
  113. y = gpuArray( single(y) );
  114. weights = gpuArray( single(weights) );
  115. for km = 1:Nvar; x0.(fitting.modelParams{km}) = gpuArray(single( x0.(fitting.modelParams{km}) ));end
  116. xStepsize = gpuArray(single(fitting.xStepSize(:)));
  117. % setup boundary variables
  118. lb = gpuArray( single(repmat(fitting.lb(:),1,Nv)));
  119. ub = gpuArray( single(repmat(fitting.ub(:),1,Nv)));
  120. % initialize array to staore all the samples
  121. xPosterior = zeros(Nvar, Nv, Ns, fitting.repetition,'single');
  123. % compute likelihood at starting points
  124. % logP is converted into external function for specific CUDA kernel
  125. % logP = @(X, Y) -sum( (this.FWD(X(1:4, :), model)-Y).^2, 1 )./(2*X(5,:).^2) + Nm/2*log(1./X(5,:).^2);
  126. xCurr = this.struct2array(x0,fitting.modelParams); % extract parameter structure to numeric array for faster computation
  127. xCurr = max(xCurr,lb); xCurr = min(xCurr,ub); % set boundary
  128. x0 = this.array2struct(xCurr,fitting.modelParams); % convert array back to structure for FWD function
  129. logP0 = arrayfun(@logP_Gaussian, sum( weights.* (FWDfunc(x0,varargin{:})-y).^2, 1 ), x0.noise, Nm);
  131. disp('-------------------------');
  132. disp('MCMC optimisation process');
  133. disp('-------------------------');
  135. % loop (multiple) proposal (same start)
  136. for ii = 1:fitting.repetition
  137. fprintf('Repetition #%i/%i \n',ii,fitting.repetition)
  139. % reset start point
  140. logPCurr = logP0;
  141. xCurr = this.struct2array(x0,fitting.modelParams);
  143. counter = 0; start = tic;
  144. for k = 1:fitting.iteration
  145. % 1. make a proposal with normal distribution
  146. % proposal is generated during iteration
  147. xProposed = xCurr + xStepsize.*randn(size(xCurr),'like',xCurr);
  148. % find proposal that is out of bound for exclusion
  149. isOutofbound = max(or(xProposed<lb, xProposed>ub),[],1);
  150. % replace out of bound by boundary values to avoid error when compting probability
  151. xProposed = max(xProposed,lb); xProposed = min(xProposed,ub);
  152. % convert the proposal into structure array for FWD function
  153. xProposed_struct = this.array2struct(xProposed,fitting.modelParams);
  155. % 2. Metropolis sampling
  156. % If the probability ratio of new to old > threshold, we take the new solution.
  157. % 2.1 proposal probability
  158. logPProposed = arrayfun(@logP_Gaussian, sum( weights.* (FWDfunc(xProposed_struct, varargin{:})-y).^2, 1 ), xProposed_struct.noise, Nm);
  159. % 2.2 Compute acceptance ratio based on new/old probability
  160. acceptanceRatio = min(exp(logPProposed-logPCurr), 1);
  161. isAccepted = acceptanceRatio > rand(1,Nv,'like',logPProposed);
  162. isAccepted(isOutofbound)= 0; % reject out of bound proposal
  163. % 2.3 update parameters if accepted
  164. logPCurr(isAccepted) = logPProposed(isAccepted);
  165. xCurr(:,isAccepted) = xProposed(:,isAccepted);
  167. % 3. Maintain the independence between iterations
  168. % 3.1 discard the first burnin*100% iterations
  169. % 3.2 keep an iteration every N iterations
  170. if ( k > Nburnin ) && mod(k-Nburnin+1, fitting.thinning) == 0 %( mod(k, fitting.thinning)==1 )
  171. counter = counter+1;
  172. xPosterior(:,:,counter,ii) = gather(xCurr);
  173. end
  175. % display message at 1000 iteration and every 10000 iteration
  176. if mod(k,fitting.iteration/50) == 0 || k == min(1e3, fitting.iteration/100)
  177. ET = duration(0,0,toc(start),'Format','hh:mm:ss');
  178. ERT = ET / (k/fitting.iteration) - ET;
  179. fprintf('Iteration #%6d, Elapsed time (hh:mm:ss):%s, Estimated remaining time (hh:mm:ss):%s \n',k,string(ET),string(ERT));
  180. end
  181. end
  182. end
  184. % convert final posterior distribution into structure
  185. xPosterior = this.array2struct(xPosterior,fitting.modelParams);
  186. for kvar = 1:Nvar; xPosterior.(fitting.modelParams{kvar}) = shiftdim(xPosterior.(fitting.modelParams{kvar}),1); end
  188. disp('The Metroplis-Hastings MCMC sampling is completed.')
  190. end
  192. function xPosterior = goodman_weare(this,y,x0,weights,fitting,modelFWD,varargin)
  193. % Input
  194. % ----------
  195. % y : measurements, [Nmeas,Nvoxels]
  196. % x0 : structure array, starting points, N fields, each field 1xNvoxel
  197. % weights : weighting for non-linear least square fitting, same dimension as y
  198. % pars0 : Structure variable containing all parameters to be estimated
  199. % fitting : Structure variable containing all fitting algorithm setting
  200. % .modelParams : 1xM cell variable, name of the model parameters, e.g. {'S0','R2star','noise'};
  201. % .lb : 1xM numeric variable, fitting lower bound, same order as field 'modelParams', e.g. [0.5, 0, 0.001];
  202. % .ub : 1xM numeric variable, fitting upper bound, same order as field 'modelParams', e.g. [2, 1, 0.1];
  203. % .iteration : # MCMC iterations
  204. % .thinning : sampling interval between iterations
  205. % .burnin : iterations at the beginning to be discarded, if burnin>1, then the exact number will be used; if 0<burnin<1 then actual burnin = iteration*burnin
  206. % .repetition : # repetition of MCMC proposal
  207. % .StepSize : step size for 'GW' in MCMC proposal ('GW' only)
  208. % .Nwalker : # random walkers ('GW' only)
  209. % .Ensembleupdate : (optional) ensemble update scheme:
  210. % 'simultaneous' - original behaviour (DEFAULT, backward compatible):
  211. % all walkers proposed/updated in one pass using a
  212. % single derangement as partners.
  213. % 'redblack' - affine-invariance-correct parallel update
  214. % (Foreman-Mackey et al. 2013): split the ensemble
  215. % into two halves and update each half using the
  216. % other (frozen) half as anchors.
  217. % FWDfunc : function handle of forward model
  218. % varargin : contains additional input requires for FWDfunc
  219. %
  220. % ENSEMBLE SAMPLERS WITH AFFINE INVARIANCE (2010) JONATHAN GOODMAN AND JONATHAN WEARE
  221. % Other references:
  222. % emcee: The MCMC Hammer (2013) https://arxiv.org/pdf/1202.3665
  223. % https://github.com/grinsted/gwmcmc/tree/master
  224. %
  226. fitting = this.check_set_default_basic(fitting);
  227. if isempty(weights); weights = ones(size(y), 'like', y); end
  229. % Nm: # measurements; Nv: # voxels
  230. [Nm, Nv] = size(y);
  231. % Nvar: # estimation parameters
  232. Nvar = numel(fitting.modelParams);
  233. Nburnin = this.get_number_burnin(fitting);
  234. % Ns: # samples in posterior distribution
  235. Ns = numel(Nburnin+1:fitting.thinning:fitting.iteration);
  236. Nwalker = fitting.Nwalker;
  237. StepSize = fitting.StepSize;
  239. % red-black update scheme
  240. useRedBlack = strcmpi(fitting.Ensembleupdate,'redblack');
  241. if useRedBlack
  242. if Nwalker < 2
  243. error('goodman_weare:Nwalker', ...
  244. 'redblack update needs Nwalker >= 2 (use even, >= 2*Nvar in practice).');
  245. end
  246. if mod(Nwalker,2) ~= 0
  247. warning('goodman_weare:Nwalker', ...
  248. 'Nwalker is odd; redblack halves will be unequal. Even Nwalker is conventional.');
  249. end
  250. h = floor(Nwalker/2);
  251. halfIdx = {1:h, h+1:Nwalker}; % {S0, S1}: complementary anchor sets
  252. end
  253. fprintf('Ensemble update scheme: %s\n', fitting.Ensembleupdate);
  255. % convert data into single datatype for better performance
  256. y = gpuArray( single(y) );
  257. weights = gpuArray( single(weights) );
  258. for km = 1:Nvar; x0.(fitting.modelParams{km}) = gpuArray(single( x0.(fitting.modelParams{km}) ));end
  260. % setup boundary variables
  261. lb = gpuArray( single(repmat(fitting.lb(:),1,Nv,Nwalker)));
  262. ub = gpuArray( single(repmat(fitting.ub(:),1,Nv,Nwalker)));
  263. % set up weight
  264. weights = repmat(weights,1,1,Nwalker);
  265. % initialize array to staore all the samples
  266. xPosterior = zeros(Nvar, Nv, Nwalker, Ns, fitting.repetition,'single');
  268. % initiate an ensemble of walkers around the starting position (0.1% full range) with Gaussian distribution
  269. % 1st: Nvar;2nd: Nv; 3rd: Nwalker
  270. xCurr = this.struct2array(x0,fitting.modelParams); % extract parameter structure to numeric array for faster computation
  271. xCurr = xCurr + (ub-lb)*fitting.startRange.*randn(size(ub)); % initiate starting position for all walkers
  272. xCurr = max(xCurr,lb); xCurr = min(xCurr,ub); % set boundary
  273. x0 = this.array2struct(xCurr,fitting.modelParams); % convert array back to structure for FWD function
  274. % compute likelihood at starting points
  275. logP0 = arrayfun(@logP_Gaussian, sum( weights.* (modelFWD(x0, varargin{:})-y).^2, 1 ), x0.noise, Nm);
  277. disp('-------------------------');
  278. disp('MCMC optimisation process');
  279. disp('-------------------------');
  281. for ii = 1:fitting.repetition
  282. fprintf('Repetition #%i/%i \n',ii,fitting.repetition)
  284. logPCurr= logP0;
  285. xCurr = this.struct2array(x0,fitting.modelParams);
  287. counter = 0; start = tic;
  288. for k = 1:fitting.iteration
  290. if ~useRedBlack
  291. % =========== LEGACY: simultaneous single-derangement update ===========
  292. % (identical to the original implementation; preserved for reproducibility)
  293. % 1.1 find a unique partner for each walker
  295. % 1. make a proposal with normal distribution
  296. % 1.1. find a unique partner for a walker k
  297. partner = this.find_partner(Nwalker);
  298. % 1.2. stretch move
  299. zz = ((StepSize-1)*rand(size(logP0),'like',xCurr) + 1).^2 / StepSize;
  300. xProposed = xCurr(:,:,partner) + (xCurr - xCurr(:,:,partner)).*zz;
  301. % find proposal that is out of bound for exclusion
  302. isOOB = max(or(xProposed<lb, xProposed>ub),[],1);
  303. % replace boundary values so it does not give error when compting probability
  304. xProposed = max(xProposed,lb); xProposed = min(xProposed,ub);
  305. % convert the proposal into structure array for FWD function
  306. xProposed_struct = this.array2struct(xProposed,fitting.modelParams);
  308. % 2. Metropolis sampling
  309. % If the probability ratio of new to old > threshold, we take the new solution.
  310. % 2.1 proposal probability
  311. logPProposed = arrayfun(@logP_Gaussian, sum( weights.* (modelFWD(xProposed_struct, varargin{:})-y).^2, 1 ), xProposed_struct.noise, Nm);
  312. % 2.2 Compute acceptance ratio based on z^(Nd-1)*new/old probability
  313. acceptanceRatio = min(zz.^(Nvar-1).*exp(logPProposed-logPCurr), 1);
  314. isAccepted = acceptanceRatio > rand(1,Nv,Nwalker,'like',xCurr);
  315. isAccepted(isOOB) = 0; % reject out of bound proposal
  316. % 2.3 update parameters
  317. logPCurr(isAccepted) = logPProposed(isAccepted);
  318. xCurr(:,isAccepted) = xProposed(:,isAccepted);
  319. else
  321. % =========== RED-BLACK: two complementary sub-steps per sweep ===========
  322. for s = 1:2
  323. active = halfIdx{s}; % walkers moved this sub-step
  324. frozen = halfIdx{3-s}; % anchors, held FIXED this sub-step
  325. na = numel(active);
  327. % stretch proposal: each active walker anchored on a RANDOM frozen
  328. % walker (uniform, WITH replacement -- textbook GW partner draw)
  329. pidx = frozen(randi(numel(frozen),1,na));
  330. xAnch = xCurr(:,:,pidx); % [Nvar,Nv,na] fixed
  331. xAct = xCurr(:,:,active); % [Nvar,Nv,na]
  332. zz = ((StepSize-1)*rand(1,Nv,na,'like',xCurr) + 1).^2 / StepSize;
  333. xProp = xAnch + (xAct - xAnch).*zz;
  335. lb_a = lb(:,:,active); ub_a = ub(:,:,active);
  336. isOOB = max(or(xProp<lb_a, xProp>ub_a),[],1);
  337. xProp = max(xProp,lb_a); xProp = min(xProp,ub_a);
  339. % proposal likelihood on the active half only
  340. % (two half-width evals/iter ~= one full eval: no cost regression)
  341. xProp_s = this.array2struct(xProp,fitting.modelParams);
  342. logPProp = arrayfun(@logP_Gaussian, ...
  343. sum( weights(:,:,active).*(modelFWD(xProp_s,varargin{:})-y).^2, 1 ), ...
  344. xProp_s.noise, Nm);
  346. % acceptance z^(Nvar-1) * pi(new)/pi(old)
  347. logPAct = logPCurr(:,:,active);
  348. aR = min(zz.^(Nvar-1).*exp(logPProp - logPAct), 1);
  349. acc = aR > rand(1,Nv,na,'like',xCurr);
  350. acc(isOOB) = 0;
  352. % scatter accepted moves back into the active block
  353. xAct(:,acc) = xProp(:,acc);
  354. logPAct(acc) = logPProp(acc);
  355. xCurr(:,:,active) = xAct;
  356. logPCurr(:,:,active) = logPAct;
  357. end
  358. end
  360. % 3. Maintain the independence between iterations
  361. % 3.1 discard the first burnin*100% iterations
  362. % 3.2 keep an iteration every N iterations
  363. if ( k > Nburnin ) && mod(k-Nburnin+1, fitting.thinning) == 0
  364. counter = counter+1;
  365. xPosterior(:,:,:,counter,ii) = gather(xCurr);
  366. end
  368. % display message at 1000 ietration and every 2000 iterations
  369. if mod(k,fitting.iteration/50) == 0 || k == min(1e2, fitting.iteration/100)
  370. ET = duration(0,0,toc(start),'Format','hh:mm:ss');
  371. ERT = ET / (k/fitting.iteration) - ET;
  372. fprintf('Iteration #%6d, Elapsed time (hh:mm:ss):%s, Estimated remaining time (hh:mm:ss):%s \n',k,string(ET),string(ERT));
  373. end
  374. end
  375. end
  377. % convert final posterior distribution into structure
  378. xPosterior = this.array2struct(xPosterior,fitting.modelParams);
  379. for kvar = 1:Nvar; xPosterior.(fitting.modelParams{kvar}) = shiftdim(xPosterior.(fitting.modelParams{kvar}),1); end
  381. disp('The affine-invariant ensemble MCMC sampling is completed.')
  383. end
  385. function xPosterior = goodman_weare_wglobal_constant(this,y,x0,weights,fitting,modelFWD,varargin)
  386. % Affine-invariant ensemble sampling for hierarchical / partial-pooling
  387. % problems with GLOBAL (shared) parameters, via block Metropolis-within-Gibbs.
  388. %
  389. % CONTRACT (scope of this sampler):
  390. % The log-likelihood must factorise over units (voxels) as
  391. % sum_z loglik(y_z | theta_z, phi)
  392. % i.e. given the global parameters phi, the units are conditionally
  393. % independent. Parameters split into LOCAL theta_z (one set per unit) and
  394. % GLOBAL phi (shared). Any model with this structure is supported
  395. % (arbitrary NND local and NGlobal global parameters); the sampler holds
  396. % no model-specific knowledge.
  397. %
  398. % WHICH PARAMETERS ARE GLOBAL is detected from the starting structure: a
  399. % field of x0 with size==1 along the voxel dimension is treated as global
  400. % (see check_global_constant). So {R2c,k2A}, a shared diffusivity, a global
  401. % B1 term, etc., all work without code changes.
  402. %
  403. % One iteration = two decoupled blocks, each a Goodman-Weare stretch move:
  404. % LOCAL block : propose theta, GLOBALS HELD FIXED, exponent z^(NND-1),
  405. % forward model evaluated at (new local, old global).
  406. % GLOBAL block: propose phi, LOCALS HELD FIXED at their updated values,
  407. % exponent z^(NGlobal-1), forward model evaluated at
  408. % (current local, new global); acceptance uses the summed
  409. % log-likelihood DIFFERENCE accumulated in double precision.
  410. %
  411. % This structure is what keeps the cached likelihood consistent with the
  412. % current state at all times (the previous single-eval / two-accept scheme
  413. % left it stale) and gives each block its correct stretch Jacobian.
  414. %
  415. % NOTES / options:
  416. % .StepSize : stretch scale 'a'
  417. % .Nwalker : # walkers (even, >= 2*max(NND,NGlobal) recommended)
  418. % .startRange : local walker init spread (fraction of range)
  419. % .globalWarmup : (optional) # global-only iterations before the joint
  420. % loop (locals frozen). Default 0. Helps when the global
  421. % posterior is very tight.
  422. % .Ensembleupdate : ensemble update scheme, applied to BOTH blocks:
  423. % 'simultaneous' - (DEFAULT, backward compatible)
  424. % all walkers proposed at once using a
  425. % single derangement as partners.
  426. % 'redblack' - correct parallel update
  427. % (Foreman-Mackey 2013): split walkers
  428. % into two halves and update each half
  429. % using the other as frozen anchors.
  430. % Consistent with the option in goodman_weare.
  431. %
  432. % Goodman & Weare 2010; emcee (Foreman-Mackey 2013).
  434. fitting = this.check_set_default_basic(fitting);
  435. if isempty(weights); weights = ones(size(y), 'like', y); end
  437. [Nm, Nv] = size(y);
  438. Nvar = numel(fitting.modelParams);
  439. Nburnin = this.get_number_burnin(fitting);
  440. Ns = numel(Nburnin+1:fitting.thinning:fitting.iteration);
  441. Nwalker = fitting.Nwalker;
  442. StepSize = fitting.StepSize;
  444. % --- ensemble update scheme (consistent with goodman_weare) -----
  445. useRedBlack = strcmpi(fitting.Ensembleupdate,'redblack');
  446. if useRedBlack
  447. if mod(Nwalker,2) ~= 0
  448. warning('goodman_weare_wglobal_constant:Nwalker', ...
  449. 'Nwalker is odd; red-black halves will be unequal. Even Nwalker is conventional.');
  450. end
  451. h = floor(Nwalker/2);
  452. halfIdx = {1:h, h+1:Nwalker};
  453. end
  454. fprintf('Ensemble update scheme: %s\n', fitting.Ensembleupdate);
  455. % -----------------------------------------------------------------
  457. % convert to single + push to GPU
  458. y = gpuArray( single(y) );
  459. weights = gpuArray( single(weights) );
  460. for km = 1:Nvar; x0.(fitting.modelParams{km}) = gpuArray(single( x0.(fitting.modelParams{km}) )); end
  462. isGlobal = this.check_global_constant(x0,fitting.modelParams);
  463. NGlobal = numel(isGlobal(isGlobal==1));
  464. NND = Nvar - NGlobal;
  465. if NND==0 || NGlobal==0
  466. error('goodman_weare_wglobal_constant:partition', ...
  467. 'Need at least one local and one global parameter (NND=%d, NGlobal=%d). Use goodman_weare for the all-local case.',NND,NGlobal);
  468. end
  470. % bounds
  471. lb = gpuArray( single(repmat(fitting.lb(isGlobal==0),1,Nv,Nwalker)));
  472. ub = gpuArray( single(repmat(fitting.ub(isGlobal==0),1,Nv,Nwalker)));
  473. lb_global = gpuArray( single(repmat(fitting.lb(isGlobal==1),1,1,Nwalker)));
  474. ub_global = gpuArray( single(repmat(fitting.ub(isGlobal==1),1,1,Nwalker)));
  475. weights = repmat(weights,1,1,Nwalker);
  477. xPosterior_ND = zeros(NND, Nv, Nwalker, Ns, fitting.repetition,'single');
  478. xPosterior_global = zeros(NGlobal, 1, Nwalker, Ns, fitting.repetition,'single');
  480. % --- initialise walkers -------------------------------------------------
  481. [xCurr_ND,xCurr_global] = this.struct2array_wglobal(x0,fitting.modelParams);
  482. % local: tight cloud around starting point
  483. xCurr_ND = xCurr_ND + (ub-lb)*fitting.startRange.*randn(size(ub));
  484. xCurr_ND = min(max(xCurr_ND,lb),ub);
  485. % GLOBAL: OVER-DISPERSED across the full prior box. Ensemble samplers
  486. % cannot manufacture spread a collapsed cloud lacks, and the global
  487. % posterior is typically very tight (informed by all voxels), so a small
  488. % cloud can get stuck. Over-dispersion contracts correctly; collapse does not.
  489. % xCurr_global = lb_global + (ub_global-lb_global).*rand(size(ub_global));
  490. xCurr_global = xCurr_global + (ub_global-lb_global)*fitting.startRange.*rand(size(ub_global));
  491. x0 = this.array2struct_wglobal(xCurr_ND,xCurr_global,fitting.modelParams,isGlobal);
  493. logP0 = arrayfun(@logP_Gaussian, sum( weights.* (modelFWD(x0,varargin{:})-y).^2, 1 ), x0.noise, Nm);
  495. disp('-------------------------');
  496. disp('MCMC optimisation process (block GW: local + global)');
  497. disp('-------------------------');
  499. for ii = 1:fitting.repetition
  500. fprintf('Repetition #%i/%i \n',ii,fitting.repetition)
  502. logPCurr = logP0;
  503. [xCurr_ND,xCurr_global] = this.struct2array_wglobal(x0,fitting.modelParams);
  505. % optional global-only warmup (locals frozen) to settle the tight global block
  506. for kw = 1:fitting.globalWarmup
  507. [xCurr_global,logPCurr] = this.gw_global_block( ...
  508. xCurr_ND,xCurr_global,logPCurr,weights,y,Nm,Nv,Nwalker, ...
  509. StepSize,NGlobal,lb_global,ub_global,fitting,isGlobal,modelFWD,varargin{:});
  510. end
  512. counter = 0; start = tic;
  513. for k = 1:fitting.iteration
  515. % ============ LOCAL block (globals fixed) ============
  516. if ~useRedBlack
  517. % --- simultaneous (default, original behaviour) ---
  518. partner = this.find_partner(Nwalker);
  519. zz_ND = ((StepSize-1)*rand(1,Nv,Nwalker,'like',xCurr_ND) + 1).^2 / StepSize;
  520. xProp_ND = xCurr_ND(:,:,partner) + (xCurr_ND - xCurr_ND(:,:,partner)).*zz_ND;
  521. oob_ND = max(or(xProp_ND<lb, xProp_ND>ub),[],1);
  522. xProp_ND = min(max(xProp_ND,lb),ub);
  523. s_local = this.array2struct_wglobal(xProp_ND, xCurr_global, fitting.modelParams, isGlobal);
  524. logP_loc = arrayfun(@logP_Gaussian, sum( weights.*(modelFWD(s_local,varargin{:})-y).^2, 1 ), s_local.noise, Nm);
  525. aR = min( zz_ND.^(NND-1).*exp(logP_loc - logPCurr), 1);
  526. acc = aR > rand(1,Nv,Nwalker,'like',xCurr_ND);
  527. acc(oob_ND) = 0;
  528. logPCurr(acc) = logP_loc(acc); % cache <- (new local, OLD global)
  529. xCurr_ND(:,acc) = xProp_ND(:,acc);
  530. else
  531. % --- red-black: two sub-steps over walker halves ---
  532. for s = 1:2
  533. active = halfIdx{s}; frozen = halfIdx{3-s}; na = numel(active);
  534. pidx = frozen(randi(numel(frozen),1,na));
  535. xAnch = xCurr_ND(:,:,pidx); xAct = xCurr_ND(:,:,active);
  536. zz_ND = ((StepSize-1)*rand(1,Nv,na,'like',xCurr_ND) + 1).^2 / StepSize;
  537. xProp = xAnch + (xAct - xAnch).*zz_ND;
  538. lb_a = lb(:,:,active); ub_a = ub(:,:,active);
  539. isOOB = max(or(xProp<lb_a, xProp>ub_a),[],1);
  540. xProp = min(max(xProp,lb_a),ub_a);
  541. % eval at (new local active, OLD global active) -- globals fixed this block
  542. s_loc = this.array2struct_wglobal(xProp, xCurr_global(:,:,active), fitting.modelParams, isGlobal);
  543. logP_loc = arrayfun(@logP_Gaussian, sum( weights(:,:,active).*(modelFWD(s_loc,varargin{:})-y).^2, 1 ), s_loc.noise, Nm);
  544. logPAct = logPCurr(:,:,active);
  545. aR = min( zz_ND.^(NND-1).*exp(logP_loc - logPAct), 1);
  546. acc = aR > rand(1,Nv,na,'like',xCurr_ND); acc(isOOB) = 0;
  547. xAct(:,acc) = xProp(:,acc);
  548. logPAct(acc) = logP_loc(acc);
  549. xCurr_ND(:,:,active) = xAct;
  550. logPCurr(:,:,active) = logPAct;
  551. end
  552. end
  554. % ============ GLOBAL block (locals fixed at updated values) ============
  555. [xCurr_global,logPCurr] = this.gw_global_block( ...
  556. xCurr_ND,xCurr_global,logPCurr,weights,y,Nm,Nv,Nwalker, ...
  557. StepSize,NGlobal,lb_global,ub_global,fitting,isGlobal,modelFWD,varargin{:});
  559. % thinning / burn-in
  560. if ( k > Nburnin ) && mod(k-Nburnin+1, fitting.thinning) == 0
  561. counter = counter+1;
  562. xPosterior_ND(:,:,:,counter,ii) = gather(xCurr_ND);
  563. xPosterior_global(:,:,:,counter,ii) = gather(xCurr_global);
  564. end
  566. if mod(k,fitting.iteration/50) == 0 || k == min(1e2, fitting.iteration/100)
  567. ET = duration(0,0,toc(start),'Format','hh:mm:ss');
  568. ERT = ET / (k/fitting.iteration) - ET;
  569. fprintf('Iteration #%6d, Elapsed time (hh:mm:ss):%s, Estimated remaining time (hh:mm:ss):%s \n',k,string(ET),string(ERT));
  570. end
  571. end
  572. end
  574. xPosterior = this.array2struct_wglobal(xPosterior_ND,xPosterior_global,fitting.modelParams,isGlobal);
  575. for kvar = 1:Nvar; xPosterior.(fitting.modelParams{kvar}) = shiftdim(xPosterior.(fitting.modelParams{kvar}),1); end
  577. disp('The block affine-invariant ensemble MCMC sampling is completed.')
  579. end
  581. % ---- GLOBAL block: one GW stretch update of phi, locals held fixed -------
  582. function [xCurr_global,logPCurr] = gw_global_block(this, ...
  583. xCurr_ND,xCurr_global,logPCurr,weights,y,Nm,Nv,Nwalker, ...
  584. StepSize,NGlobal,lb_global,ub_global,fitting,isGlobal,modelFWD,varargin)
  586. useRedBlack = strcmpi(fitting.Ensembleupdate,'redblack');
  587. if ~useRedBlack
  588. % --- simultaneous (default, original behaviour) ---
  589. partner = this.find_partner(Nwalker);
  590. zz_g = ((StepSize-1)*rand(1,1,Nwalker,'like',xCurr_global) + 1).^2 / StepSize;
  591. gProp = xCurr_global(:,:,partner) + (xCurr_global - xCurr_global(:,:,partner)).*zz_g;
  592. oob_g = max(or(gProp<lb_global, gProp>ub_global),[],1);
  593. gProp = min(max(gProp,lb_global),ub_global);
  594. s_glob = this.array2struct_wglobal(xCurr_ND, gProp, fitting.modelParams, isGlobal);
  595. logP_g = arrayfun(@logP_Gaussian, sum( weights.*(modelFWD(s_glob,varargin{:})-y).^2, 1 ), s_glob.noise, Nm);
  596. dlogP = sum( double(logP_g - logPCurr), 2 );
  597. aR_g = min( zz_g.^(NGlobal-1).*exp(dlogP), 1 );
  598. acc_g = aR_g > rand(1,1,Nwalker,'like',xCurr_global);
  599. acc_g(oob_g) = 0; acc_g = logical(acc_g);
  600. xCurr_global(:,:,acc_g) = gProp(:,:,acc_g);
  601. logPCurr(:,:,acc_g) = logP_g(:,:,acc_g);
  602. else
  603. % --- red-black: two sub-steps over walker halves ---
  604. h_g = floor(Nwalker/2); halfIdx_g = {1:h_g, h_g+1:Nwalker};
  605. for s = 1:2
  606. active = halfIdx_g{s}; frozen = halfIdx_g{3-s}; na = numel(active);
  607. pidx = frozen(randi(numel(frozen),1,na));
  608. xg_act = xCurr_global(:,:,active);
  609. zz_g = ((StepSize-1)*rand(1,1,na,'like',xCurr_global) + 1).^2 / StepSize;
  610. gProp = xCurr_global(:,:,pidx) + (xg_act - xCurr_global(:,:,pidx)).*zz_g;
  611. oob_g = max(or(gProp<lb_global(:,:,active), gProp>ub_global(:,:,active)),[],1);
  612. gProp = min(max(gProp,lb_global(:,:,active)),ub_global(:,:,active));
  613. % eval at (current local active, new global active) -- locals fixed this block
  614. s_glob = this.array2struct_wglobal(xCurr_ND(:,:,active), gProp, fitting.modelParams, isGlobal);
  615. logP_g = arrayfun(@logP_Gaussian, sum( weights(:,:,active).*(modelFWD(s_glob,varargin{:})-y).^2, 1 ), s_glob.noise, Nm);
  616. lp_act = logPCurr(:,:,active);
  617. dlogP = sum( double(logP_g - lp_act), 2 ); % [1,1,na], double
  618. aR_g = min( zz_g.^(NGlobal-1).*exp(dlogP), 1 );
  619. acc_g = aR_g > rand(1,1,na,'like',xCurr_global);
  620. acc_g(oob_g) = 0; acc_g = logical(acc_g);
  621. xg_act(:,:,acc_g) = gProp(:,:,acc_g);
  622. lp_act(:,:,acc_g) = logP_g(:,:,acc_g);
  623. xCurr_global(:,:,active) = xg_act;
  624. logPCurr(:,:,active) = lp_act;
  625. end
  626. end
  627. end
  629. end
  631. methods(Static)
  633. % check and set default fitting algorithm parameters
  634. function fitting2 = check_set_default_basic(fitting)
  635. % Input
  636. % -----
  637. % fitting : structure contains fitting algorithm parameters
  638. % .iteration : no. of maximum MCMC iterations, default = 200k
  639. % .repetition : no. of MCMC repetitions, default = 1
  640. % .thinning : MCMC thinning interval, default = every 20 iterations
  641. % .burnin : MCMC burn-in ratio, default = 10%
  642. % .metric : method to compute expected valur from posterior distribution, 'mean' (default) | 'median'
  643. % .algorithm : MCMC algorithm 'MH': Metropolis-Hastings; 'GW': Goodman-Weare, 'MH' (default) | 'GW'
  644. % .StepSize : Step size for Goodman-Weare, default = 2
  645. % .Nwalker : number of walkers for Goodman-Weare, default = 50
  646. %
  647. fitting2 = fitting;
  649. % get fitting algorithm setting
  650. if ~isfield(fitting,'iteration'); fitting2.iteration = 2e5; end
  651. if ~isfield(fitting,'thinning'); fitting2.thinning = 20; end % thinning, sampled every 100 interval
  652. if ~isfield(fitting,'metric'); fitting2.metric = {'mean','std'}; end
  653. if ~isfield(fitting,'burnin'); fitting2.burnin = 0.1; end % 10% burnin
  655. if ~isfield(fitting,'repetition'); fitting2.repetition = 1; end
  656. if ~isfield(fitting,'outputFilename'); fitting2.outputFilename = []; end
  657. if ~isfield(fitting,'algorithm'); fitting2.algorithm = 'MH'; end
  658. if ~isfield(fitting,'StepSize'); fitting2.StepSize = 2; end
  659. if ~isfield(fitting,'Nwalker'); fitting2.Nwalker = 50; end
  660. if ~isfield(fitting,'ub'); fitting2.ub = []; end
  661. if ~isfield(fitting,'lb'); fitting2.lb = []; end
  662. if ~isfield(fitting,'startRange'); fitting2.startRange = 0.001; end
  664. % --- choose ensemble update scheme (default = original behaviour) ----
  665. if ~isfield(fitting,'Ensembleupdate'); fitting2.Ensembleupdate = 'simultaneous'; end
  666. if ~isfield(fitting,'globalWarmup'); fitting2.globalWarmup = 0; end
  668. if any(ismember(fitting2.metric,'mode'))
  669. if ~isfield(fitting,'Nbin'); fitting2.Nbin = 1001; end
  670. end
  672. if ~isfield(fitting,'autoMemManage'); fitting2.autoMemManage = true; end
  675. if ~iscell(fitting2.metric)
  676. fitting2.metric = cellstr(fitting2.metric);
  677. end
  679. if strcmpi(fitting2.algorithm ,'gw'); fitting2.algorithm = 'ensemble'; end % legacy
  680. end
  682. % display fitting algorithm parameters
  683. function display_basic_algorithm_parameters(fitting)
  685. if strcmpi( fitting.algorithm, 'ensemble'); algorithm = 'Affine-Invariant Ensemble';
  686. else; algorithm = 'Metropolis-Hastings'; end
  688. disp('----------------------------------------------------');
  689. disp('Markov Chain Monte Carlo (MCMC) algorithm parameters');
  690. disp('----------------------------------------------------');
  691. disp(['Algorithm : ', algorithm]);
  692. disp(['No. of iterations : ', num2str(fitting.iteration)]);
  693. disp(['No. of repetitions: ', num2str(fitting.repetition)])
  694. disp(['Thinning : ', num2str(fitting.thinning)]);
  695. disp(['Burn-in (#iter.) : ' num2str(mcmc.get_number_burnin(fitting))])
  696. disp(['Metric(s) : ', cell2str(fitting.metric)]);
  697. if strcmpi( fitting.algorithm, 'ensemble'); disp(['Step size : ', num2str(fitting.StepSize) ]); end
  698. if strcmpi( fitting.algorithm, 'ensemble'); disp(['No. of walkers : ', num2str(fitting.Nwalker) ]); end
  700. end
  702. % save the mcmc output structure variable into disk space
  703. function save_mcmc_output(outputFilename,out)
  704. % Input
  705. % ------------------
  706. % outputFilename : output filename
  707. % out : output structure of askadam
  708. %
  710. % save the estimation results if the output filename is provided
  711. if ~isempty(outputFilename)
  712. [output_dir,~,~] = fileparts(outputFilename);
  713. if ~exist(output_dir,'dir')
  714. mkdir(output_dir);
  715. end
  716. save(outputFilename,'out');
  717. fprintf('Estimation output is saved at %s\n',outputFilename);
  718. end
  719. end
  721. % convert numerical array into structure variable for FWD function
  722. function x_struct = array2struct(x,fields)
  723. for k = 1:numel(fields)
  724. x_struct.(fields{k}) = x(k,:,:,:,:,:);
  725. end
  726. end
  728. % convert structure variable into numerical array
  729. function x = struct2array(x_struct,fields)
  731. nVol = size(x_struct.(fields{1}),2);
  732. nWalker = size(x_struct.(fields{1}),3);
  734. x = gpuArray(zeros(numel(fields),nVol,nWalker,"single"));
  735. for k = 1:numel(fields)
  736. x(k,:,:) = x_struct.(fields{k});
  737. end
  738. end
  740. function x_struct = array2struct_wglobal(x_ND,x_global,fields, isGlobal)
  742. ctr_global = 1; ctr_ND = 1;
  743. for k = 1:numel(fields)
  744. if isGlobal(k)
  745. x_struct.(fields{k}) = x_global(ctr_global,:,:,:,:,:);
  746. ctr_global = ctr_global + 1;
  747. else
  748. x_struct.(fields{k}) = x_ND(ctr_ND,:,:,:,:,:);
  749. ctr_ND = ctr_ND +1;
  750. end
  751. end
  752. end
  754. function [x_ND,x_global] = struct2array_wglobal(x_struct,fields)
  756. % check if the fitting parameter is global or voxel
  757. isGlobal = mcmc.check_global_constant(x_struct,fields);
  758. nVol = 0; for k = 1:numel(fields); nVol = max(size(x_struct.(fields{k}),2),nVol); end
  760. nWalker = size(x_struct.(fields{1}),3);
  761. x_ND = gpuArray(zeros(numel(isGlobal(isGlobal==0)),nVol,nWalker,"single"));
  762. x_global = gpuArray(zeros(numel(isGlobal(isGlobal==1)),1,nWalker,"single"));
  764. ctr_global = 1; ctr_ND = 1;
  765. for k = 1:numel(fields)
  766. if isGlobal(k)
  767. x_global(ctr_global,1,:) = x_struct.(fields{k});
  768. ctr_global = ctr_global +1;
  769. else
  770. x_ND(ctr_ND,:,:) = x_struct.(fields{k});
  771. ctr_ND = ctr_ND +1;
  772. end
  773. end
  774. end
  776. function isGlobal = check_global_constant(x_struct,fields)
  777. % check if the fitting parameter is global or voxel
  778. isGlobal = zeros(numel(fields),1);
  779. for k = 1:numel(fields)
  780. if size(x_struct.(fields{k}),2) > 1
  781. isGlobal(k) = false;
  782. else
  783. isGlobal(k) = true;
  784. end
  785. end
  786. end
  788. % compute the number of iteration requires for burn-in
  789. function Nburnin = get_number_burnin(fitting)
  790. if fitting.burnin < 1
  791. Nburnin = floor(fitting.iteration*fitting.burnin);
  792. else
  793. Nburnin = fitting.burnin;
  794. end
  795. end
  797. % find a unique partner for an index
  798. function partner = find_partner(maxIndex)
  799. isSelfPartner = true;
  800. while isSelfPartner
  801. partner = randperm(maxIndex);
  802. isSelfPartner = any(partner == 1:maxIndex,'all');
  803. end
  804. end
  806. % convert estimation into organised output structure
  807. function out = res2out(xPosterior,fitting,mask)
  809. % store the unshaped posterior into out
  810. out.posterior = xPosterior;
  812. % compute additional metric if specified
  813. fields = fieldnames(xPosterior);
  815. % Nvox = size(xPosterior.(fields{1}),1);
  816. Nsample = prod(size(xPosterior.(fields{1}),2:5));
  818. metrics = fitting.metric;
  819. if ~isempty(metrics)
  820. for km = 1:numel(metrics)
  821. switch lower(metrics{km})
  822. case 'mean'
  823. for kvar=1:numel(fields)
  824. tmp = mean( reshape( xPosterior.(fields{kvar}), [size(xPosterior.(fields{kvar}),1), Nsample]),2);
  825. tmp = utils.reshape_ND2image(tmp,mask);
  826. out.mean.(fields{kvar}) = tmp;
  827. end
  828. case 'median'
  829. for kvar=1:numel(fields)
  830. tmp = median( reshape( xPosterior.(fields{kvar}), [size(xPosterior.(fields{kvar}),1), Nsample]),2);
  831. tmp = utils.reshape_ND2image(tmp,mask);
  832. out.median.(fields{kvar}) = tmp;
  833. end
  834. case 'std'
  835. for kvar=1:numel(fields)
  836. tmp = std( reshape( xPosterior.(fields{kvar}), [size(xPosterior.(fields{kvar}),1), Nsample]),[],2);
  837. tmp = utils.reshape_ND2image(tmp,mask);
  838. out.std.(fields{kvar}) = tmp;
  839. end
  840. case 'iqr'
  841. for kvar=1:numel(fields)
  842. tmp = iqr( reshape( xPosterior.(fields{kvar}), [size(xPosterior.(fields{kvar}),1), Nsample]),2);
  843. tmp = utils.reshape_ND2image(tmp,mask);
  844. out.iqr.(fields{kvar}) = tmp;
  845. end
  846. case 'mode'
  847. for kvar=1:numel(fields)
  849. Nbin = fitting.Nbin;
  851. idx = find(ismember(fitting.modelParams,fields{kvar}));
  852. edges = linspace(fitting.lb(idx)-1e-8,fitting.ub(idx)+1e-8,Nbin);
  854. tmp = reshape( xPosterior.(fields{kvar}), [size(xPosterior.(fields{kvar}),1), Nsample]);
  855. tmp = mode(discretize(tmp,edges),2);
  856. tmp = (edges(tmp) + edges(tmp+1)) / 2;
  857. tmp = utils.reshape_ND2image(tmp.',mask);
  858. out.mode.(fields{kvar}) = tmp;
  859. end
  860. end
  861. end
  862. end
  863. end
  865. % make sure all network parameters stay between 0 and 1
  866. function parameters = set_boundary(parameters,ub,lb)
  868. field = fieldnames(parameters);
  869. for k = 1:numel(field)
  870. parameters.(field{k}) = max(parameters.(field{k}),lb(k)); % Lower bound
  871. parameters.(field{k}) = min(parameters.(field{k}),ub(k)); % upper bound
  873. end
  875. end
  877. end
  878. end

mcmc.m at commit 05f49d4, under GPL-3.0 · at the source

Overview

Authors: Yixin Ma1,2, Laleh Eskandarian1,2, Kwok‐Shing Chan1,2, Hansol Lee1,2,3, Gabriel Ramos‐Llordén1,2, Aneri Bhatt1, Julianna Gerold1, Mirsad Mahmutovic4, Boris Keil4,5, Hong‐Hsi Lee1,2, Susie Y Huang1,2,6
ORCID iDs: Kwok‐Shing Chan, Hansol Lee, Aneri Bhatt, Mirsad Mahmutovic, Boris Keil, Susie Y Huang
  1. Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Charlestown, Massachusetts, USA
  2. Harvard Medical School, Boston, Massachusetts, USA
  3. Department of Biomedical Engineering, Ulsan National Institute of Science and Technology, Ulsan, South Korea
  4. Institute of Medical Physics and Radiation Protection, Mittelhessen University of Applied Sciences, Giessen, Germany
  5. Department of Diagnostic and Interventional Radiology, University Hospital Marburg, Philipps University of Marburg, Marburg, Germany
  6. Harvard‐MIT Division of Health Sciences and Technology, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
Institutions: Harvard University (United States); Massachusetts General Hospital (United States); Athinoula A. Martinos Center for Biomedical Imaging (United States); Ulsan National Institute of Science and Technology (South Korea); Technische Hochschule Mittelhessen (Germany); Philipps University of Marburg (Germany); Universitätsklinikum Gießen und Marburg (Germany); Harvard–MIT Division of Health Sciences and Technology (United States); Massachusetts Institute of Technology (United States)
Journal: Human brain mapping, volume 47, issue 8, article e70553
Dates: received 29 August 2025; accepted 30 April 2026; published online 4 June 2026; in print June 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1002/hbm.70553 · PMID 42240067 · PMCID PMC13266422 · OpenAlex W7163552186
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism)
Methods: Spectral & time-frequency, Statistics, Connectivity, fMRI & imaging
MeSH: Axons*, Brain*, Connectome*, Diffusion Magnetic Resonance Imaging*, Adult, Female, Humans, Male, Middle Aged, Young Adult (* major topic)
Topic: Advanced Neuroimaging Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: National Institute on Aging (R21AG085795, K99AG073506); NIBIB NIH HHS (U01EB026996-S1, P41EB030, P41EB015896, U01EB026996); National Institute of Biomedical Imaging and Bioengineering (U01EB026996‐S1, U01EB026996, P41EB030, P41EB015896); NINDS NIH HHS (U24NS137077, R01NS118187); Office of the Director (OD) of the National Institutes of Health (S10OD032184); ZonMw Rubicon Fellowship (04520232330012); Office of the Director (OD) of the NIH in partnership with the National Institute of Dental & Craniofacial Research (NIDCR) (DP5OD031854); NIA NIH HHS (K99AG073506, R21AG085795); National Institute of Neurological Disorders and Stroke (U24NS137077, R01NS118187); ZonMw (04520232330012); NIH/NINDS Pathway-to-Independence Award (K99NS132984); National Research Foundation of Korea (NRF) (RS-2024-00411768)
Citations: not cited yet (Europe PMC); 89 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 5 matches between paragraphs and lines of code.

yixinma9/diffusion_preproc_C2

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 59032d4f5ea70d3fbe37d93ae8d4d049e217d315, 7 August 2025
Languages: Python (15), Shell (1)
Size: 29 files, 16 scripts
Software Heritage: not archived
Found in: the text, “Image Preprocessing”
Holds: README, environment (requirements.txt)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NiBabel (4 files), MRtrix3 (3 files), NumPy (3 files), FSL (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
17 files
At the source: github.com/yixinma9/diffusion_preproc_C2

kschan0214/gacelle

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 05f49d48dac7926b6c69a5612ec5dfca7f8cee1b, 9 September 2026
Languages: MATLAB (210), Python (1)
Size: 276 files, 211 scripts
Software Heritage: not archived
Found in: the text, “Markov Chain Monte Carlo With Ensemble Sampler W”
Holds: README, license file, environment (docs/requirements.txt), tests, continuous integration, documentation
Not found: CITATION.cff
Tools: Parallel Computing Toolbox (30 files), Image Processing Toolbox (28 files), Deep Learning Toolbox (12 files), Optimization Toolbox (10 files), FreeSurfer (1 file), Statistics and Machine Learning Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
213 files
At the source: github.com/kschan0214/gacelle

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 227 scripts, each with its path and the digest of its content;
  • 5 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1002/hbm.70553.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 10 MeSH terms, 12 funders, 87 references.

Cite

This paper

Ma, Y., Eskandarian, L., Chan, K., Lee, H., Ramos‐Llordén, G., Bhatt, A., Gerold, J., Mahmutovic, M., Keil, B., Lee, H., & Huang, S. Y. (2026). Axon Diameter Mapping in the Living Human Brain with Ultra-High-Gradient Diffusion MRI at 500 mT/m Gradient Strength. Human brain mapping, 47(8), e70553. https://doi.org/10.1002/hbm.70553

BibTeX

@article{ma2026axon,
author = {Ma, Yixin and Eskandarian, Laleh and Chan, Kwok‐Shing and Lee, Hansol and Ramos‐Llordén, Gabriel and Bhatt, Aneri and Gerold, Julianna and Mahmutovic, Mirsad and Keil, Boris and Lee, Hong‐Hsi and Huang, Susie Y},
title = {{Axon Diameter Mapping in the Living Human Brain with Ultra-High-Gradient Diffusion MRI at 500 mT/m Gradient Strength}},
journal = {Human brain mapping},
year = {2026},
month = jun,
volume = {47},
number = {8},
pages = {e70553},
publisher = {Wiley},
issn = {1065-9471},
doi = {10.1002/hbm.70553},
url = {https://doi.org/10.1002/hbm.70553},
pmid = {42240067},
pmcid = {PMC13266422}
}

RIS

TY - JOUR
AU - Ma, Yixin
AU - Eskandarian, Laleh
AU - Chan, Kwok‐Shing
AU - Lee, Hansol
AU - Ramos‐Llordén, Gabriel
AU - Bhatt, Aneri
AU - Gerold, Julianna
AU - Mahmutovic, Mirsad
AU - Keil, Boris
AU - Lee, Hong‐Hsi
AU - Huang, Susie Y
TI - Axon Diameter Mapping in the Living Human Brain with Ultra-High-Gradient Diffusion MRI at 500 mT/m Gradient Strength
T2 - Human brain mapping
J2 - Hum Brain Mapp
PY - 2026
DA - 2026/06/01
VL - 47
IS - 8
SP - e70553
SN - 1065-9471
PB - Wiley
DO - 10.1002/hbm.70553
UR - https://doi.org/10.1002/hbm.70553
LA - en
ER -

CSL-JSON

{
"id": "10.1002/hbm.70553",
"type": "article-journal",
"title": "Axon Diameter Mapping in the Living Human Brain with Ultra-High-Gradient Diffusion MRI at 500 mT/m Gradient Strength",
"container-title": "Human brain mapping",
"author": [
{
"family": "Ma",
"given": "Yixin"
},
{
"family": "Eskandarian",
"given": "Laleh"
},
{
"family": "Chan",
"given": "Kwok‐Shing"
},
{
"family": "Lee",
"given": "Hansol"
},
{
"family": "Ramos‐Llordén",
"given": "Gabriel"
},
{
"family": "Bhatt",
"given": "Aneri"
},
{
"family": "Gerold",
"given": "Julianna"
},
{
"family": "Mahmutovic",
"given": "Mirsad"
},
{
"family": "Keil",
"given": "Boris"
},
{
"family": "Lee",
"given": "Hong‐Hsi"
},
{
"family": "Huang",
"given": "Susie Y"
}
],
"container-title-short": "Hum Brain Mapp",
"volume": "47",
"issue": "8",
"page": "e70553",
"DOI": "10.1002/hbm.70553",
"PMID": "42240067",
"PMCID": "PMC13266422",
"ISSN": "1065-9471",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/hbm.70553",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1109/tmi.2026.3664328 [code]
Axon Diameter Mapping From Myelin Water Diffusion MRI.
Journal: IEEE transactions on medical imaging
In common: MRtrix3, Optimization Toolbox, Image Processing Toolbox, 1 other tool, structural MRI / diffusion, 23 references, 2 authors
[2] doi:10.1002/mrm.70490 [code]
Dependence of the Extra-Cellular Diffusion Coefficient on the Fractions of Neurites and Cell Bodies in Gray Matter.
Journal: Magnetic resonance in medicine
In common: structural MRI / diffusion, 17 references, author Hansol Lee
[3] doi:10.1002/mrm.70378 [code]
Investigating the Sensitivity of the Diffusion MRI Signal to Magnetization Transfer and Permeability via Monte-Carlo Simulations.
Journal: Magnetic resonance in medicine
In common: FSL, NumPy, structural MRI / diffusion, 16 references
[4] doi:10.1371/journal.pbio.3003861
Learning engages transient and sustained cellular mechanisms in the human brain.
Journal: PLoS biology
In common: structural MRI / diffusion, 12 references
[5] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: MRtrix3, Optimization Toolbox, Parallel Computing Toolbox, 6 other tools, 1 reference
[6] doi:10.1162/imag.a.1317 [code]
Bayesian insights into exchange and restriction in gray matter diffusion MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: NiBabel, NumPy, structural MRI / diffusion, 8 references
[7] doi:10.1038/s41598-026-51531-w [code]
Multimodal age-dependent diffusion-MRI analysis of the neocortex in a rat model of cortical dysplasia.
Journal: Scientific reports
In common: MRtrix3, FSL, Image Processing Toolbox, 3 other tools, structural MRI / diffusion, 4 references
[8] doi:10.7554/elife.107661 [code]
In vivo mapping of striatal neurodegeneration in Huntington's disease with Soma and Neurite Density Imaging.
Journal: eLife
In common: Deep Learning Toolbox, Optimization Toolbox, Image Processing Toolbox, 1 other tool, structural MRI / diffusion, 4 references
[9] doi:10.1038/s43856-026-01614-6 [code]
Simulation-based inference at the theoretical limit for fast, robust microstructural MRI with minimal diffusion data.
Journal: Communications medicine
In common: NumPy, structural MRI / diffusion, 7 references
[10] doi:10.1038/s41586-026-10631-3 [code]
A prognostic human brain network for diffuse midline glioma.
Journal: Nature
In common: Optimization Toolbox, Parallel Computing Toolbox, FreeSurfer, 5 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.