-
Notifications
You must be signed in to change notification settings - Fork 8
/
nip_spm_priors.m
105 lines (102 loc) · 3.72 KB
/
nip_spm_priors.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
function [h,Qp] = nip_spm_priors(y,L,patches,Qe,type)
% function h = nip_spm_priors(y,L,patches,Qe,type)
% Find the weights for a set of cortical "patches" using ARD or GS as
% implemented in the SPM software package
% Input:
% y -> NcxNt. Matrix containing the data,
% L -> NcxNd. Lead Field matrix
% Patches -> NdxNp. Each column of the matrix contains the each of
% the spatial pattern.
% Qe -> NcxNc. Covariance matrix at the sensor level.
% type -> string. Method used to find the hyperparameters:
% 'ARD' - 'GS'
% Output:
% h -> Npx1. Weights for each of the patches.
%
% Additional Comments: This function is based on a script written by
% Jose David Lopez - [email protected]
% Gareth Barnes - [email protected]
% Vladimir Litvak - [email protected]
% This function is a wrapper to some SPM utilities!
%
% Juan S. Castaño
% 9. May 2013
% YY = y*y';
[Nd, Np] = size(patches);
% Nc = size(L,1);
% Qp = {};
% LQpL = {};
%
% eyeNd = repmat(eye(Nd),[1,1,3]);
%
% Lnew = zeros(Nc,Nc);
Qp = zeros(3*Nd,Np*3);
for i = 1:Np
% i
if size(L,2)~= size(patches,1);
Ltemp = nip_translf(L);
for k = 1:3
for j = 1:3
if k==j
% Lnew(:,:) = Ltemp(:,:,j)*diag(patches(:,i))*Ltemp(:,:,j)';
Q(:,1,j) = patches(:,i);
else
Q(:,1,j) = 0.0*ones(Nd,1,1);
end
end
% LQpL{end+1}.q = Lnew;
Qp(:,(i-1)*3+k) = nip_translf(permute(Q,[2 1 3]))';
% Qp{end+1}.q = nip_translf(permute(Q,[2 1 3]))';
Q = zeros(Nd,1,3);
end
else
% Qp{end + 1}.q = patches(:,i);
% LQpL{end + 1}.q = L*diag(patches(:,i))*L';
end
end
% Np = numel(Qp);
switch(type)
case {'GS'}
% Greedy search over MSPs
% Np = length(Qp);
% Q = zeros(3*Nd,numel(Qp));
% for i = 1:size(Q,2)
% Q(:,i) = Qp{i}.q;
% end
% Q = sparse(Q);
% Q = patches;
% Multivariate Bayes (Here is performed the inversion)
%------------------------------------------------------------------
MVB = spm_mvb(y,L,[],Qp,Qe,16);
% Accumulate empirical priors (New set of patches for the second inversion)
%------------------------------------------------------------------
% MVB.cp provides the final weights of the hyperparameters
% Qcp = Q*MVB.cp;
% QP{end + 1} = sum(Qcp.*Q,2);
% LQP{end + 1} = (L*Qcp)*Q';
% LQPL{end + 1} = LQP{end}*L';
h = diag(MVB.cp);
case {'ARD'}
% ReML - ARD (Here is performed the inversion)
%------------------------------------------------------------------
[Cy,h,Ph,F] = spm_sp_reml(YY,[],[Qe LQpL],1);
h = h(2:end);
% Spatial priors (QP)
%------------------------------------------------------------------
% h provides the final weights of the hyperparameters
% Ne = length(Qe);
% Np = length(Qp);
% hp = h(Ne + (1:Np));
% qp = sparse(0);
% for i = 1:Np
% if hp(i) > max(hp)/128;
% qp = qp + hp(i)*Qp{i}.q*Qp{i}.q';
% end
% end
%
% % Accumulate empirical priors (New set of patches for the second inversion)
% %------------------------------------------------------------------
% QP{end + 1} = diag(qp);
% LQP{end + 1} = L*qp;
% LQPL{end + 1} = LQP{end}*L';
end