OSCR

Joint estimation of source dynamics and interactions from MEG data

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] § MATERIALS AND METHODS › Reference Source Estimation With Beamformers and MVAR Fitting ↔ joint_est/saem.m, lines 1–34 · score 0.65 · lead field matrix, covariance matrix, source space, iterative, Ns, MCMV
  2. [2] § MATERIALS AND METHODS › Real MEG Data › Preprocessing. ↔ EcoG_simu/extract_ecog_data.m, lines 24–102 · score 0.63 · pass filtered, preprocessed, Epochs, baseline, fT, 0.3 Hz
  3. [3] § MATERIALS AND METHODS › Simulation Using ECoG Signals ↔ joint_est/saem.m, lines 1–34 · score 0.60 · source space, measurement noise, MVAR Models, vertices, algorithms, Ns
  4. [4] § MATERIALS AND METHODS › Reference Source Estimation With Beamformers and MVAR Fitting ↔ two_step_est/source_localizers/mkfilt_lcmv.m, the whole file · a weak match · score 0.60 · lead field matrix, covariance matrix, LCMV, beamformers, filters, Ns
  5. [5] § MATERIALS AND METHODS › JEDI-MEG ↔ joint_est/saem.m, lines 419–499 · score 0.50 · Kalman filter, source amplitudes, PMCMC, particle, SAEM

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 · 501 lines · 17 KB · MIT · 3 matches

  1. function [state, params, LL_complete]= saem(init, model, opt_params, Y_avg, Y, Y_mcmv, Ns, Np, P, dist_thr, gpu_flag, whiten_flag, mcmv_prior_flag)
  2. %AUTHOR : NARAYAN SUBRAMANIYAM /AALTO/NBE/
  3. % INPUTS
  4. % 1) init : struct containing initial values for parameters and kalman
  5. % mean/cov . init.A0, init.V_q0, init.q0, init.P0
  7. % 2) model : struct containing model params
  8. % model.G : lead-field matrix
  10. %model.XI : measurement noise covariance matrix
  12. %model.V_r : jitter matrix for particles
  14. %model.mesh : mesh details (mesh.p, mesh.e)
  16. %model.incl_vert : list of vertices included in source space. Its dimension
  17. %should be the same as the no. of columns in lead field
  21. %model.anat_parc : anatomical parcels (this will be optional in future.
  22. %Presently this info is used to initialize the particles in different
  23. %anatomical areas
  25. % 3) Y : MEG data (channels X T)
  27. % 4) opt_params : struct containing parameters the algorithm depends on
  28. %opt_params.Ns : no of sources / dipoles
  29. %opt_params.p : MVAR model order
  30. %opt_params.N : number of particles
  31. %opt_params.num_iter : no of iterations
  32. %opt_params.gamma : forgetting factor for SAEM
  33. %opt_params.dist_thr : distance threshold (just to make sure initial particles
  34. %draw are not too close to each other. This can be some value like 2 or 3
  37. %%%%%%%%%%%%%%PULL OUT VALUES FROM STRUCTS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  40. % source space, lead-field matrix
  41. G = model.G;
  42. pnts_e = model.pnts_e;
  43. incl_vert = model.incl_verts; % source points used to calculate lead field
  44. pnts = pnts_e(incl_vert, :); %xyz of vertices included in source-space
  47. %initial parameter values
  48. A_tilde = init.A0; %AR matrix
  49. V_tilde = init.V0;
  50. P_tilde = init.P0;
  51. q_tilde = init.q0;
  53. if whiten_flag==0
  54. sigma_m = init.sigma_m0;
  55. sigma_b = init.sigma_b0;
  56. %measurement+bio noise cov matrix
  57. E = sigma_m^2*eye(size(G,1)) + sigma_b^2*(G*G');
  58. else
  59. E = eye(size(G,1));
  60. end
  64. % optimization parameters
  65. Niter = opt_params.Niter;
  66. gamma = opt_params.gamma;
  69. % set few constants
  70. Nq = Ns*P;
  71. Ny = size(Y,1);
  72. T = size(Y,2);
  73. Nj=size(Y,3); %no of trials
  74. Nr=Ns*3;
  75. CONST = -Ny * 0.5 * log(2*pi);
  77. % initialize sufficient statistics matrix
  78. S = zeros(3*Nq+Ny, 3*Nq+Ny);
  84. % parameter and state histories
  85. params.A_hist = zeros(Ns, Ns*P, Niter+1);
  86. params.A_hist(:,:,1) = A_tilde(1:Ns,:);
  87. params.V_hist(:,:,1) = V_tilde(1:Ns, 1:Ns);
  88. if whiten_flag==0
  89. params.sigma_m_hist(1) = sigma_m;
  90. params.sigma_b_hist(1) = sigma_b;
  91. end
  92. state.r_hat_est = zeros(Nr, T);
  93. state.qs_hat_est = zeros(Ns, T+1, Nj);
  94. state.qf_hat_est = zeros(Ns, T, Nj);
  95. state.r_hat_hist = zeros(Nr,Niter+1);
  99. %log-likelihoods
  100. LL_complete = zeros(Niter,1);
  102. %some variables for gpu
  103. if gpu_flag ==1
  104. look_up = zeros(Nj, Np);
  105. for j=1:Nj
  106. look_up(j,:) = (j-1)*Np+1 : Np*j;
  107. end
  108. Y_rep = zeros(Ny, Np, Nj);
  109. end
  113. % Random initialization of conditional trajectory
  115. if mcmv_prior_flag==0
  116. [ind_prime, r_prime, G_prime] = set_part_traj...
  117. (pnts, G, dist_thr, Ns, T); %for bbcb_simu
  118. else
  119. [ind_prime, r_prime, G_prime] = set_part_traj_mcmv_prior(Y_mcmv, eye(size(Y_mcmv,1)), G, pnts_e, incl_vert, Ns, size(Y_mcmv,2));
  121. end
  122. state.r_hat_hist(:,1) = mean(r_prime,2);
  123. m0=mean(r_prime,2);
  124. P0 = 16*eye(3*Ns);
  127. for k=1:1:Niter
  128. sprintf('iteration no %d', k)
  129. tic
  131. %Kalman filter variables. We assume Kalman covariance is same across
  132. %trials to keep things computationally feasible.
  134. q_tilde_kf = single(zeros(Nq,T,Np,Nj));
  135. q_tilde_kp = single(zeros(Nq,T,Np,Nj));
  136. P_tilde_kf = single(zeros(Nq,Nq,T,Np));
  137. P_tilde_kp = single(zeros(Nq,Nq,T,Np));
  139. % PMCMC variables
  140. r_pf = zeros(Nr,T,Np);
  141. ind_pf = zeros(Ns, T, Np);
  142. G_pf=zeros(Ny,Ns,T,Np);
  144. % weights
  145. w = single(zeros(Np,T));
  146. log_weights = single(zeros(Np,T));
  148. % replace Np-th particle with conditional trajectory
  149. r_pf(:,:,end) = r_prime(:,1:T);
  150. G_pf(:,:,:,end) = G_prime(:,:,1:T);
  151. ind_pf(:,:,end) = ind_prime(:,1:T);
  153. % initialize Kalman filter
  154. q_tilde_kp(:,1,:,:) = repmat(q_tilde,[1,1,Np, Nj]);
  155. P_tilde_kp(:,:,1,:) = repmat(P_tilde,[1,1,1, Np]);
  157. %initialize weights
  158. w(:,1) = ones(size(w(:,1)));
  159. w(:,1) = w(:,1)./sum(w(:,1));
  161. % stats for backward smoothing
  162. omega_tr = zeros(Ns, Ns, T-1, Nj);
  163. lambda_tr = zeros(Ns, T-1, Nj);
  166. % Ancestor indices
  167. an_hist = zeros(Np, T);
  170. % Compute lambda, omega for conditional trajectory
  172. for j=1:1:Nj
  173. [omega_tr(:,:,:,j), lambda_tr(:,:,j)] = compute_omega_lambda_fast...
  174. (A_tilde(1:Ns,1:Ns), V_tilde, G_pf, squeeze(Y(:,:,j)), E, Ns, P, Np, T);
  175. end
  177. omega_avg=squeeze(mean(omega_tr,4));
  178. lambda_avg=squeeze(mean(lambda_tr,3));
  180. % draw particles for t=1
  181. [ind_pf(:,1,1:end-1), r_pf(:,1,1:end-1), G_pf(:,:,1,1:end-1)] = ...
  182. set_all_part_traj_v2(pnts, G, dist_thr, Ns, Np-1,m0,P0);
  183. % draw particles for t=1
  184. %[ind_pf(:,1,1:end-1), r_pf(:,1,1:end-1), G_pf(:,:,1,1:end-1)] = ...
  185. %set_all_part_traj(pnts, G, dist_thr, Ns, Np-1);
  188. % stochastic E-STEP
  189. disp('entering time loop')
  190. for t=1:T
  192. if t >=2
  193. [a,~] = find(mnrnd(1,repmat(w(:,t-1)',Np-1,1))');
  194. [ind_pf(:,t,1:end-1), XI] = LW_model(squeeze(r_pf(:,t-1,a)), 0.95, pnts);
  195. log_det = 2*sum(log(diag(chol(XI))));
  197. for i=1:1:Ns
  198. r_pf(3*(i-1)+1:3*i,t,1:end-1) = pnts(ind_pf(i,t,1:end-1),:)';
  199. end
  201. for i=1:1:Np-1
  202. G_pf(:,:,t,i) = G(:,squeeze(ind_pf(:,t,i)));
  203. end
  204. q_tilde_kf_mean = squeeze(mean(q_tilde_kf(1:Ns,t-1,:,:),4));
  205. P_tilde_kf_mean = squeeze(P_tilde_kf(1:Ns,1:Ns,t-1,:))./Nj;
  206. if gpu_flag == 1
  207. parallel.gpu.enableCUDAForwardCompatibility(1)
  208. [aN] = compute_aN_gpu(repmat(r_pf(:,t,end), [1 1 Np]), r_pf(:,t-1,:), q_tilde_kf_mean, P_tilde_kf_mean, XI, omega_avg(1:Ns,1:Ns,t-1), lambda_avg(1:Ns,t-1),...
  209. w(:,t-1), Np, CONST, log_det);
  210. else
  211. [aN] = compute_aN_fast(repmat(r_pf(:,t,end), [1 1 Np]), r_pf(:,t-1,:), q_tilde_kf_mean, P_tilde_kf_mean, XI, omega_avg(:,:,t-1), lambda_avg(:,t-1),...
  212. Ns, w(:,t-1), Np, CONST, log_det);
  213. end
  214. an_hist(:,t) = [a;aN];
  215. end
  217. if gpu_flag == 1 % perform kalman predict, update and log-likelihood computation on GPU
  218. G_tilde = squeeze(G_pf(:,:,t,:));
  219. for j=1:Nj
  220. Y_rep(:,:,j) = repmat(squeeze(Y(:,t,j)), [1, 1, Np]);
  221. end
  222. Y_reshape = reshape(Y_rep, Ny, Np*Nj);
  224. %kalman predict(GPU)
  225. if t >= 2
  226. % tic
  227. q_tilde_kf_ = reshape(squeeze(q_tilde_kf (:,t-1,:,:)), Nq, Np*Nj);
  228. [q_tilde_kp_, P_tilde_kp(:,:,t,:)] = gpu_kalman_predict(q_tilde_kf_, squeeze(P_tilde_kf(:,:,t-1,:)), A_tilde, V_tilde);
  229. q_tilde_kp (:,t,:,:) = reshape(q_tilde_kp_, Nq, Np, Nj);
  230. % toc
  231. end
  233. % kalman update(GPU)
  234. %tic
  235. q_tilde_kp_ = reshape(squeeze(q_tilde_kp(:,t,:,:)), Nq, Np*Nj);
  236. q_tilde_kp_gpu = gpuArray(q_tilde_kp_);
  237. P_tilde_kp_gpu = gpuArray(squeeze(P_tilde_kp(:,:,t,:)));
  238. Y_gpu = gpuArray(Y_reshape);
  240. G_tilde_ = repmat(G_tilde, [1 1 1 Nj]);
  241. G_tilde_=reshape(G_tilde_, Ny,Ns,Np*Nj);
  242. G_tilde_gpu = gpuArray(G_tilde);
  243. G_tilde_gpu_ = gpuArray(G_tilde_);
  244. E_gpu = gpuArray(E);
  245. [q_tilde_kf_, P_tilde_kf(:,:,t,:), Kf_] = gpu_kalman_update(q_tilde_kp_gpu, P_tilde_kp_gpu, Y_gpu, G_tilde_gpu, G_tilde_gpu_, E_gpu, Np, Nj);
  246. q_tilde_kf(:,t,:,:) = reshape(q_tilde_kf_, Nq, Np, Nj);
  247. %toc
  248. for i=1:Np
  249. tmp_G_tilde = squeeze(G_tilde(:,:,i));
  250. ypred(:,i) = tmp_G_tilde*mean(q_tilde_kp_(1:Ns, look_up(:,i)),2);
  251. sigma(:,:,i) = (tmp_G_tilde*squeeze(P_tilde_kp (1:Ns,1:Ns,t,i))*tmp_G_tilde' + E )./Nj;
  252. log_det_sigma(i) = 2*sum(log(diag(chol(squeeze(sigma(:,:,i))))));
  253. end
  254. %tic
  255. log_weights(:,t) = gpu_logpdf(Y_avg(:,t),ypred,sigma, log_det_sigma,Np,Ny);
  256. %toc
  257. else
  258. for i = 1:Np
  259. for j=1:1:Nj
  260. % Prediction
  261. if t >= 2
  262. [q_tilde_kp(:,t,i,j), P_tilde_kp(:,:,t,i)] = kalman_predict(q_tilde_kf(:,t-1,i,j), P_tilde_kf(:,:,t-1,i), A_tilde, V_tilde );
  264. end
  266. % Update
  267. tmp_G_tilde = [G_pf(:,:,t,i) repmat(zeros(size(G_pf(:,:,t,i))), 1, P-1)];
  268. [q_tilde_kf(:,t,i,j), P_tilde_kf(:,:,t,i), Kf(:,:,i)] = kalman_update(q_tilde_kp(:,t,i,j), P_tilde_kp(:,:,t,i), squeeze(Y(:,t,j)),tmp_G_tilde, E);
  270. end
  271. % log-likelihood
  272. ypred(:,i) = tmp_G_tilde*squeeze(mean(q_tilde_kp(:,t,i,:),4));
  273. sigma(:,:,i) = (tmp_G_tilde*P_tilde_kp(:,:,t,i)*tmp_G_tilde' + E )./Nj;
  274. log_weights(i,t)=logpdf(Y_avg(:,t), ypred(:,i), sigma(:,:,i));
  275. end
  276. end
  278. % PF weight update
  279. maxlog = max(log_weights(:,t));
  280. log_weights(:,t) = log_weights(:,t) - maxlog;
  281. w(:,t) = exp(log_weights(:,t));
  282. w(:,t) = w(:,t) / sum(w(:,t));
  283. end
  285. % set trajectories based on ancestral history
  286. ind_an = an_hist(:,T);
  287. for t = T-1:-1:1
  288. r_pf(:,t,:) = r_pf(:,t,ind_an);
  289. q_tilde_kp(:,t,:,:) = q_tilde_kp(:,t,ind_an,:);
  290. P_tilde_kp(:,:,t,:) = P_tilde_kp(:,:,t,ind_an);
  291. q_tilde_kf(:,t,:,:) = q_tilde_kf(:,t,ind_an,:);
  292. P_tilde_kf(:,:,t,:) = P_tilde_kf(:,:,t,ind_an);
  293. G_pf(:,:,t,:) = G_pf(:,:,t,ind_an);
  294. ind_an = an_hist(ind_an,t);
  295. end
  300. if gpu_flag==1
  301. q_tilde_ks = zeros(Nq,T+1,Np*Nj); P_tilde_ks = zeros(Nq,Nq,T+1,Np);
  302. M = zeros(Nq,Nq,T,Np);
  304. q_tilde_kp = reshape(q_tilde_kp, Nq, T, Np*Nj);
  305. q_tilde_kf = reshape(q_tilde_kf, Nq, T, Np*Nj);
  307. P_tilde_ks(:,:,end,:) = squeeze(P_tilde_kp(:,:,end,:)); P_tilde_ks(:,:,end-1,:) = squeeze(P_tilde_kf(:,:,end,:));
  308. q_tilde_ks(:,end,:) = squeeze(q_tilde_kp(:,end,:)); q_tilde_ks(:,end-1,:) = squeeze(q_tilde_kf(:,end,:));
  310. P_tilde_ks_t1 = squeeze(P_tilde_ks(:,:,end-1,:));
  312. q_tilde_ks_t1 = squeeze( q_tilde_ks(:,end-1,:));
  316. for t = T-1:-1:1
  317. P_tilde_kf_t = squeeze(P_tilde_kf(:,:,t,:));
  318. P_tilde_kp_t1 = squeeze(P_tilde_kp(:,:,t+1,:));
  319. q_tilde_kf_t = squeeze(q_tilde_kf(:,t,:));
  320. q_tilde_kp_t1 = squeeze(q_tilde_kp(:,t+1,:));
  321. [J(:,:,t,:), q_tilde_ks(:,t,:), P_tilde_ks(:,:,t,:)] ...
  322. = gpu_kalman_smoother(P_tilde_kf_t, P_tilde_kp_t1, q_tilde_kf_t,q_tilde_kp_t1, ...
  323. q_tilde_ks_t1, P_tilde_ks_t1, A_tilde, Np, Nj);
  324. P_tilde_ks_t1 = squeeze(P_tilde_ks(:,:,t,:));
  325. q_tilde_ks_t1 = squeeze( q_tilde_ks(:,t,:));
  326. end
  330. % Computing of P_{t-1,t|T}
  331. G_T = [G_pf(:,:,T,:) repmat(zeros(size(G_pf(:,:,T,:))), 1, P-1)];
  332. G_T=squeeze(G_T);
  333. I=eye(Nq);
  334. M(:,:,end,:) = gpu_init_M(I, squeeze(Kf_(:,:,:)), G_T, squeeze(P_tilde_kf(:,:,T-1,:)), A_tilde);
  335. for t=T-1:-1:2
  336. P_tilde_kf_t = squeeze(P_tilde_kf(:,:,t,:));
  337. M_tp1 = squeeze(M(:,:,t+1,:));
  338. J_tm1 = squeeze(J(:,:,t-1,:));
  339. J_t = squeeze(J(:,:,t,:));
  340. M(:,:,t,:) = gpu_one_lag(A_tilde, P_tilde_kf_t, J_tm1, M_tp1, J_t);
  341. end
  342. q_tilde_ks = reshape(q_tilde_ks, Nq, T+1, Np, Nj);
  343. q_tilde_kf = reshape(q_tilde_kf, Nq, T, Np, Nj);
  344. else
  345. q_tilde_ks = zeros(Nq,T+1,Np,Nj); P_tilde_ks = zeros(Nq,Nq,T+1,Np);
  346. M = zeros(Nq,Nq,T,Np);
  348. for i = 1:Np
  349. % Inititalizing
  350. P_tilde_ks(:,:,end,i) = P_tilde_kp(:,:,end,i); P_tilde_ks(:,:,end-1,i) = P_tilde_kf(:,:,end,i);
  351. for t = T-1:-1:1
  352. J(:,:,t,i) = P_tilde_kf(:,:,t,i)*A_tilde'/P_tilde_kp(:,:,t+1,i);
  353. P_tilde_ks(:,:,t,i) = P_tilde_kf(:,:,t,i) + J(:,:,t,i)*(P_tilde_ks(:,:,t+1,i)-P_tilde_kp(:,:,t+1,i))*J(:,:,t,i)';
  354. end
  355. for j=1:Nj
  356. q_tilde_ks(:,end,i,j) = q_tilde_kp(:,end,i,j); q_tilde_ks(:,end-1,i,j) = q_tilde_kf(:,end,i,j);
  357. for t = T-1:-1:1
  358. q_tilde_ks(:,t,i,j) = q_tilde_kf(:,t,i,j)+J(:,:,t,i)*(q_tilde_ks(:,t+1,i,j)-q_tilde_kp(:,t+1,i,j));
  359. end
  360. end
  361. % Computing of M = P_{t-1,t|T}
  362. tmp_ = [G_pf(:,:,T,i) repmat(zeros(size(G_pf(:,:,T,i))), 1, P-1)];
  363. M(:,:,end,i) = (eye(Nq)-Kf(:,:,i)*tmp_)*A_tilde*P_tilde_kf(:,:,T-1,i);
  365. for t = T-1:-1:2
  366. M(:,:,t,i) = P_tilde_kf(:,:,t,i)*J(:,:,t-1,i)' + J(:,:,t,i)*(M(:,:,t+1,i)-A_tilde*P_tilde_kf(:,:,t,i))*J(:,:,t-1,i)';
  367. end
  368. end
  370. end
  372. star = logical(mnrnd(1,w(:,end)));
  373. r_prime = squeeze(r_pf(:,:,star));
  374. ind_prime = squeeze(ind_pf(:,:,star));
  375. G_prime = squeeze(G_pf(:,:,:,star));
  378. % M-STEP
  379. %Bt_s = zeros(3*Nq+Ny,3*Nq+Ny, T-1);
  382. % compute sufficient stats
  384. S3T = zeros(3*Nq+Ny,3*Nq+Ny);
  387. % compute sufficient stats
  388. Y_tmp = repmat(Y,[1 1 1 Np]);
  390. for t=2:T
  391. w_ = permute(repmat(w(:,t), [1 3*Nq+Ny 3*Nq+Ny 1]), [2 3 1]);
  392. w_1= squeeze(w_(:,1,:));
  393. Bt = compute_Bt(P_tilde_ks(:,:,t,:), P_tilde_ks(:,:,t-1,:), M(:,:,t,:), Nq, Ny, Np);
  394. Bt_s = sum(w_.*Bt, 3);
  395. for j=1:1:Nj
  396. S3T = suff_stats_fast(S3T, q_tilde_ks(:,t,:,j), q_tilde_ks(:,t-1,:,j), squeeze(Y_tmp(:,t,j,:)), Bt_s, Np, Nq, Ny, w_1);
  397. end
  398. end
  400. S = (1-gamma(k)).*S + gamma(k).*S3T;
  401. S=double(S);
  402. Phi = S(1:Nq,1:Nq);
  403. Psi = S(1:Nq,Nq+(1:Nq));
  404. Sigma = S(Nq+(1:Nq),Nq+(1:Nq));
  405. Xi_1 = S(2*Nq+Ny+(1:Nq),2*Nq+Ny+(1:Nq));
  406. Lambda_1 = S(2*Nq+(1:Ny),2*Nq+Ny+(1:Nq));
  407. Omega_1 = S(2*Nq+(1:Ny),2*Nq+(1:Ny));
  408. G_s = mean(squeeze(mean(G_pf, 4)),3);
  412. stats.Phi = Phi(1:Ns, 1:Ns);
  413. stats.Psi = Psi(1:Ns, 1:Nq);
  414. stats.Sigma = Sigma;
  415. stats.Xi = Xi_1(1:Ns,1:Ns);
  416. stats.Lambda = Lambda_1(:,1:Ns);
  417. stats.Omega = Omega_1;
  419. %compute complete likelihood before optimization
  420. LLT = compute_complete_ll(Nj*T, stats, V_tilde(1:Ns, 1:Ns), A_tilde(1:Ns,1:Nq), E, G_s(:,1:Ns));
  422. %estimate V_q
  423. V_tilde = (Phi - (Psi/Sigma)*Psi')/(Nj*T);
  424. V_tilde(Ns+1:Ns*P, 1:Ns)=0;
  425. V_tilde(1:Ns*P, Ns+1:end)=0;
  426. V_tilde(1:Ns, 1:Ns)=V_tilde(1:Ns, 1:Ns).*eye(Ns);
  428. % check if V_q is SPD
  429. V_tilde(1:Ns, 1:Ns) = max(V_tilde(1:Ns, 1:Ns),1e-4);
  430. if min(eig(V_tilde(1:Ns, 1:Ns)))< 0
  431. V_tilde(1:Ns, 1:Ns) = V_tilde(1:Ns, 1:Ns) + 0.001*eye(Ns);
  432. end
  433. V_tilde(1:Ns, 1:Ns);
  434. %estimate source amplitudes
  435. % q_ks_hat = mean(squeeze(mean(q_tilde_ks,4)),3);
  436. % q_kf_hat = mean(squeeze(mean(q_tilde_kf,4)),3);
  438. % estimate A
  439. tmp = Psi/Sigma;
  440. A_tilde = tmp(:,1:Nq);
  441. A_tilde(Ns+1:Ns*P, 1:Ns*(P-1)) = eye(Ns*(P-1));
  442. A_tilde(Ns+1:Ns*P, 1+Ns*(P-1):end) = zeros(Ns*(P-1), Ns);
  445. %estimate E
  446. if whiten_flag == 0
  447. S_R = Omega_1 - Lambda_1(:,1:Ns)*G_s' - G_s*Lambda_1(:,1:Ns)' + G_s*Xi_1(1:Ns,1:Ns)*G_s';
  448. sigma_m = grad_descent_sigmam(S_R, G, Nj*T, sigma_b, sigma_m)
  449. sigma_b = grad_descent_sigmab(S_R, G, Nj*T, sigma_b, sigma_m)
  450. E = sigma_m^2*eye(size(G,1)) + sigma_b^2*(G*G');
  451. % check if E is SPD
  452. [V_,D_] = eig(E);
  453. [I_,J_] = find(D_<0);
  455. if (~isempty(I_))
  456. for ii = 1:length(I_)
  457. D_(I_(ii),J_(ii)) = 1e-11;
  458. end
  459. E = V_*D_*V_';
  460. end
  463. end
  465. %compute complete likelihood after optimization
  466. LLTopt = compute_complete_ll(Nj*T, stats, V_tilde(1:Ns, 1:Ns), A_tilde(1:Ns,1:Nq), E, G_s(:,1:Ns));
  468. %store the histories of estimated parameters
  469. params.A_hist(:,:,k+1) = A_tilde(1:Ns,:);
  470. params.V_hist(:,:,k+1) = V_tilde(1:Ns, 1:Ns);
  471. if whiten_flag==0
  472. params.sigma_b_hist(1,k+1) = sigma_b;
  473. params.sigma_m_hist(1,k+1) = sigma_m;
  474. end
  475. display(['Iteration ',num2str(k),'. Increase in LL: ', num2str(LLTopt - LLT)])
  476. LL_complete(k,1) = LLTopt;
  479. % compute posterior means from PMCMC and Kalman filter
  480. r_pf_hat = zeros(Nr, T);
  481. for t=1:1:T
  482. r_pf_hat(:,t) = sum(repmat(w(:,t),...
  483. 1,3*Ns) .* ...
  484. squeeze(r_pf(:,t,:))');
  485. end
  489. state.r_hat_est = r_pf_hat;
  490. state.r_hat_hist(:,k) = mean(r_pf_hat,2);
  491. state.qf_hat_est = squeeze(mean(q_tilde_kf(1:Ns,:,:,:),3));
  492. state.qs_hat_est = squeeze(mean(q_tilde_ks(1:Ns,:,:,:),3));
  493. state.r_particles = r_pf;
  494. %state.qs_hist(:,:) = q_ks_hat(1:Ns,:);
  495. %state.qs_trials = squeeze(mean(q_tilde_ks(1:Ns,:,:,:),3));
  496. %state.qf_hist(:,:) = q_kf_hat(1:Ns,:);
  497. toc
  498. end
  499. end

saem.m at commit 3e0189c, under MIT · at the source

Overview

Authors: Narayan Puthanmadam Subramaniyam1, Filip Tronarp2, Simo Särkkä2, Lauri Parkkonen1
ORCID iDs: Lauri Parkkonen
  1. Department of Neuroscience and Biomedical Engineering, Aalto University, Espoo, Finland
  2. Department of Electrical Engineering and Automation, Aalto University, Espoo, Finland
Institutions: Aalto University (Finland)
Journal: n/a, volume 9, issue 3, pages 842-868
Dates: received 5 September 2024; accepted 4 March 2025; published online 17 July 2025
Type: Research article · Language: English
License: none stated
Identifiers: DOI 10.1162/netn_a_00453 · PMCID PMC12283153 · OpenAlex W4408559139
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: MEG (modality)
Methods: Spectral & time-frequency, Smoothing, state filtering, decompositions, Preprocessing, Statistics, Physiology & signal measures
Keywords: MEG, Bayesian filtering, Functional connectivity, Source localization
Topic: Underwater Acoustics Research (Oceanography, Earth and Planetary Sciences), according to OpenAlex
Funding: Research Council of Finland (289108)
Citations: not cited yet (Europe PMC); 45 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (none stated) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

narayanps/jediMEG

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 3e0189cc3eb18a4613c875bfdcf872fe644f12dd, 26 January 2025
Languages: MATLAB (102)
Size: 127 files, 102 scripts
Software Heritage: not archived
Found in: “SUPPORTING INFORMATION”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Parallel Computing Toolbox (9 files), Statistics and Machine Learning Toolbox (3 files), FieldTrip (1 file), Signal Processing Toolbox (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
104 files
At the source: github.com/narayanps/jediMEG

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 102 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

Datasets cited

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, 4 authors, 4 keywords, 1 funder, 43 references.

Cite

This paper

Puthanmadam Subramaniyam, N., Tronarp, F., Särkkä, S., & Parkkonen, L. (2025). Joint estimation of source dynamics and interactions from MEG data. Network Neuroscience, 9(3), 842-868. https://doi.org/10.1162/netn_a_00453

BibTeX

@article{puthanmadamsubramaniyam2025joint,
author = {Puthanmadam Subramaniyam, Narayan and Tronarp, Filip and Särkkä, Simo and Parkkonen, Lauri},
title = {{Joint estimation of source dynamics and interactions from MEG data}},
journal = {Network Neuroscience},
year = {2025},
volume = {9},
number = {3},
pages = {842--868},
publisher = {MIT Press},
issn = {2472-1751},
doi = {10.1162/netn_a_00453},
url = {https://doi.org/10.1162/netn_a_00453},
pmcid = {PMC12283153}
}

RIS

TY - JOUR
AU - Puthanmadam Subramaniyam, Narayan
AU - Tronarp, Filip
AU - Särkkä, Simo
AU - Parkkonen, Lauri
TI - Joint estimation of source dynamics and interactions from MEG data
T2 - Network Neuroscience
J2 - Netw Neurosci
PY - 2025
DA - 2025
VL - 9
IS - 3
SP - 842
EP - 868
SN - 2472-1751
PB - MIT Press
DO - 10.1162/netn_a_00453
UR - https://doi.org/10.1162/netn_a_00453
LA - en
ER -

CSL-JSON

{
"id": "10.1162/netn_a_00453",
"type": "article-journal",
"title": "Joint estimation of source dynamics and interactions from MEG data",
"container-title": "Network Neuroscience",
"author": [
{
"family": "Puthanmadam Subramaniyam",
"given": "Narayan"
},
{
"family": "Tronarp",
"given": "Filip"
},
{
"family": "Särkkä",
"given": "Simo"
},
{
"family": "Parkkonen",
"given": "Lauri"
}
],
"container-title-short": "Netw Neurosci",
"volume": "9",
"issue": "3",
"page": "842-868",
"DOI": "10.1162/netn_a_00453",
"PMCID": "PMC12283153",
"ISSN": "2472-1751",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/netn_a_00453",
"language": "en",
"issued": {
"date-parts": [
[
2025
]
]
}
}

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.1002/hbm.70484 [code]
Frequency-Resolved Cortical Functional Connectivity Across the Adult Lifespan.
Journal: Human brain mapping
In common: MEG, 3 references, author Lauri Parkkonen
[2] doi:10.1038/s41597-025-06397-4 [code]
MEG-SCANS - A comprehensive magnetoencephalography speech dataset with Stories, Chirps and Noisy Sentences
Journal: n/a
In common: FieldTrip, Parallel Computing Toolbox, Signal Processing Toolbox, 1 other tool, MEG, 1 reference
[3] doi:10.1038/s41467-026-71604-8 [code]
Respiration as a dynamic modulator of sensory sampling.
Journal: Nature communications
In common: FieldTrip, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, MEG, 2 references
[4] doi:10.1007/s12021-025-09759-w [code]
Limitations of Variational Laplace-Based Dynamic Causal Modelling for Multistable Cortical Circuits.
Journal: Neuroinformatics
In common: Parallel Computing Toolbox, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 2 references
[5] doi:10.1002/advs.77857 [code]
Brain Network Dynamics of Local and Global Predictive Processing in Aging.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: FieldTrip, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 2 references
[6] doi:10.1111/nyas.70349 [code]
FREQ-NESS Reveals Age-Related Differences in Frequency-Resolved Brain Networks During Auditory Recognition and Resting State.
Journal: Annals of the New York Academy of Sciences
In common: FieldTrip, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, MEG, 1 reference
[7] doi:10.1002/hbm.70624 [code]
Exploring the Sensitivity Limits of Neuronal Current Imaging With MRI and MEG in the Human Brain.
Journal: Human brain mapping
In common: FieldTrip, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, MEG, 1 reference
[8] doi:10.1162/imag.a.1219 [code]
Causal mechanisms of individual differences in hemispheric lateralization of the face perception network: A DCM-PEB approach.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: 4 references
[9] doi:10.1371/journal.pone.0353371 [code]
Inter-brain functional connectivity: Are we measuring the right thing?
Journal: PloS one
In common: MEG, author Lauri Parkkonen
[10] doi:10.1038/s41593-026-02205-3 [code]
Competitive interactions shape mammalian brain network dynamics and computation.
Journal: Nature neuroscience
In common: Parallel Computing Toolbox, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 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.