Skip to content

Commit 1781a19

Browse files
author
prni2017-smc
committed
Added source code
1 parent 8ae4d31 commit 1781a19

18 files changed

+1293
-0
lines changed

‎src/+covariance/LICENSE.md

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
Copyright (C) 2017 Manjari Narayan
2+
3+
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
4+
5+
1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
6+
7+
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8+
9+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

‎src/+covariance/check_symposdef.m

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
function [issym isposdef results] = check_symposdef(A)
2+
3+
[p1 p2] = size(A);
4+
assert(p1==p2,'Not square matrix');
5+
6+
vecA = A(find(triu(ones(p1,p2),1)));
7+
vecB = A(find(tril(ones(p1,p2),-1)));
8+
sym_err = sum((vecA-vecB).^2)/length(vecA);
9+
issym = sym_err<.5;
10+
if(issym & sym_err> 1e-2)
11+
warning('Numerical Asymmetry. Symmetricize by A + A^T');
12+
issym = 1;
13+
else
14+
issym = 2;
15+
end
16+
17+
[V D] = eig(A);
18+
tol = 1e-5;
19+
if(min(D)>tol)
20+
isposdef = 2; % Positive Definite
21+
elseif(abs(min(D))<=tol)
22+
isposdef = 1; % Positive Semi-Definite
23+
else
24+
isposdef = 0; % Negative Definite
25+
end
26+
27+
results.sym_err = sym_err;
28+
results.eigs = D;
29+
results.eigv = V;
30+
end

‎src/+covariance/mle_sample_covariance.m

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
function [SigHat results] = mle_sample_covariance(X,varargin)
2+
%MLE_SAMPLE_COVARIANCE
3+
%
4+
5+
6+
if(ndims(X)==3)
7+
[m p n] = size(X);
8+
elseif(ndims(X)==2)
9+
[m p] = size(X);
10+
n = 1;
11+
end
12+
13+
narginchk(1,2);
14+
if(nargin==2)
15+
options = varargin{1};
16+
if(~isfield(options,'M'))
17+
options.M = eye(m);
18+
end
19+
if(~isfield(options,'verbose'))
20+
options.verbose = true;
21+
end
22+
if(~isfield(options,'standardize'))
23+
options.standardize = true;
24+
end
25+
else
26+
% default
27+
options.M = eye(m);
28+
options.verbose = true;
29+
options.standardize = true;
30+
end
31+
32+
if(options.standardize);
33+
[Y succnorm] = standardize.successive_normalize(X');
34+
X = Y';
35+
else
36+
[X succnorm] = standardize.standardize_cols(X);
37+
end
38+
39+
% Initialize
40+
SigHat = zeros(p,p);
41+
M = options.M;
42+
43+
if(M(1,1)~=1)
44+
M = M/M(1,1);
45+
if(options.verbose)
46+
sprintf('Scaling Weight Matrix by first entry')
47+
end
48+
end
49+
50+
if n>1
51+
for cc=1:n
52+
SigHat = SigHat + X(:,:,cc)'*M*X(:,:,cc)/m;
53+
end
54+
SigHat = SigHat/n;
55+
else
56+
SigHat = (X'*M*X)/m;
57+
end
58+
59+
results.X = X;
60+
results.whitenMatrix = M;
61+
results.succnorm = succnorm;
62+
results.options = options;
63+
64+
end

‎src/+covariance/sample_conditional_correlation.m

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
1+
function [Sigma results] = sample_conditional_correlation(X,options)
2+
%SAMPLE_CONDITIONAL_CORRELATION
3+
%
4+
% Sigma = Sigma_obs - Sigma_nuisance
5+
%
6+
% Obtained by
7+
% Schur complement using covariance decomposition
8+
% Residuals after regressing out confounds and modified cholesky decomposition
9+
%
10+
% USAGE:
11+
%
12+
% INPUT
13+
% - X : A single n x p Data matrix
14+
% - options.scaled: Columns of X scaled to unit diagonals
15+
% - options.verbose: Plot and print outputs
16+
% - options.outputdir: Defaults to 'tmp/'
17+
% - options.filename: Defaults to 'conditionalCov'
18+
%
19+
% OUTPUT
20+
% - results.Sigma - Denoised Correlation Matrix
21+
% - results.nCov - Nuisance Covariance Matrix
22+
% - results.NSR - Nuisance to Observed Signal Ratio (trace(nCov)/trace(Sigma_obs))
23+
% - results.X_perpY - Projection of Data onto orthogonal complement of nuisance Y
24+
% - results.nCorr - Rescaled Nuisance Correlation Matrix
25+
26+
if(~exist('options','var'))
27+
options = []
28+
end
29+
if(isempty(options))
30+
options = []
31+
end
32+
if(~isfield(options,'verbose'))
33+
options.verbose = false;
34+
end
35+
if(~isfield(options,'outputdir'))
36+
options.outputdir = fullfile('tmp',datestr(now,'dd-mmm-yyyy-HHMM'));
37+
end
38+
if(~isfield(options,'filename'))
39+
options.filename = mfilename;
40+
end
41+
42+
% TBD: Add process_options or processArgs
43+
% ~isfield(options,'verbose')
44+
45+
if(~isfield(options,'corrfun'))
46+
%options.corrfun = @(X)(cov(X));
47+
options.corrfun = @(X)(covariance.mle_sample_covariance( ...
48+
X, ...
49+
struct('standardize',false) ...
50+
));
51+
end
52+
53+
function Sig = my_xcorr(X,y,options)
54+
% Returns p x q cross correlation between n x p and n x q
55+
56+
[n1 p] = size(X);
57+
[n2 q] = size(y);
58+
assert(n1==n2,'X and y need to have same row length')
59+
assert(n2>q,'Too many nuisance variables');
60+
61+
if(~exist('options','var'))
62+
options = [];
63+
end
64+
if(~isfield(options,'standardize'))
65+
options.standardize = false;
66+
end
67+
68+
if(options.standardize)
69+
% Xt = standardize.successive_normalize(X');
70+
% X = Xt';
71+
X = standardize.standardize_cols(X);
72+
if(q>1)
73+
%y = standardize.successive_normalize(y);
74+
[~, mu_y, sig_y] = standardize.standardize_cols(y);
75+
else
76+
[~, mu_y, sig_y] = standardize.standardize_cols(y);
77+
end
78+
Sig = X'*y/n1;
79+
sig_y(sig_y<=0) = 1;
80+
Sig = Sig*diag(1./sig_y);
81+
82+
else
83+
[~, mu_x, sig_x] = standardize.standardize_cols(X);
84+
[~, ~, sig_y] = standardize.standardize_cols(y);
85+
%y = zscore(y);
86+
87+
if(q>1)
88+
Sig = bsxfun(@minus,X,mu_x)'*y/n1;
89+
elseif(q==1)
90+
Sig = bsxfun(@minus,X,mu_x)'*y/n1;
91+
end
92+
93+
sig_x(sig_x<=0) = 1;
94+
Sig = diag(1./sig_x)*Sig;
95+
%Sig = diag(1./sig_x)*Sig*diag(1./sig_y);
96+
97+
end
98+
99+
end
100+
101+
proj = @(X)(X*pinv(X'*X)*X');
102+
if(~isfield(options,'xcorrfun'))
103+
%options.xcorrfun = @(X,y)(my_xcorr(X,y,struct('standardize',true)));
104+
options.xcorrfun = @(X,y)(proj(y)*standardize.standardize_cols(X));
105+
end
106+
107+
if(~isfield(options,'nuisance'))
108+
options.nuisance = mean(X,2);
109+
end
110+
111+
y = options.nuisance;
112+
Sig_xx = options.corrfun(X);
113+
Sig_xy = options.xcorrfun(X,y);
114+
Sig_yy = options.corrfun(y);
115+
if(size(Sig_xy,2)==size(Sig_yy,1))
116+
tmpSigmaL = Sig_xy*pinv(Sig_yy)*Sig_xy';
117+
[LV LD] = eig(tmpSigmaL);
118+
[sortLD, LDidx] = sort(diag(LD),'descend'); LDidx = LDidx(real(sortLD)>0);
119+
SigmaL = LV(:,LDidx)*real(LD(LDidx,LDidx))*LV(:,LDidx)';
120+
SigmaL = real(SigmaL);
121+
else
122+
tmpSigmaL = cov(Sig_xy);
123+
SigmaL = tmpSigmaL;
124+
end
125+
126+
% Remove low rank component
127+
Sigma = Sig_xx - SigmaL;
128+
varcorr = @(D,Sig)(diag(1./D)*Sig*diag(1./D));
129+
D = real(diag(Sigma));
130+
D(D<=0) = 1;
131+
132+
results.SigXX = Sig_xx;
133+
results.SigXY = Sig_xy;
134+
results.SigYY = Sig_yy;
135+
results.nCov = SigmaL;
136+
results.Sigma = varcorr(sqrt(D),Sigma);
137+
results.SigmaCov = Sigma; %
138+
results.nCorr = varcorr(sqrt(diag(SigmaL)),SigmaL);
139+
results.NSR = trace(SigmaL)/trace(Sig_xx);
140+
results.SNR = trace(Sigma)/trace(SigmaL);
141+
results.Y = y;
142+
143+
function Xy = Xy_orthogonalize(X,Y)
144+
% Nuisance signal regression
145+
146+
% p = size(X,2);
147+
% Xy = zeros(size(X));
148+
%
149+
% for ii=1:p
150+
% gsProj = gsr'*X(:,ii)/(gsr'*gsr)*gsr;
151+
% Xy = X(:,ii) - gsProj;
152+
% end
153+
154+
mu = mean(X);
155+
Xcen = bsxfun(@minus,X,mu);
156+
157+
proj = @(X)(X*pinv(X'*X)*X');
158+
Xy = Xcen - proj(Y)*Xcen;
159+
Xy = bsxfun(@plus,Xy,mu);
160+
161+
end
162+
163+
results.X_perpY = Xy_orthogonalize(standardize.standardize_cols(X),options.nuisance);
164+
165+
166+
if(options.verbose)
167+
% Only perform for comparison purposes;
168+
nuisance = options.nuisance;
169+
Xgsr = Xy_orthogonalize(standardize.standardize_cols(X),nuisance);
170+
171+
if(exist('brewermap'))
172+
colormapfun = @()(brewermap(length(colormap),'RdYlBu'));
173+
close all;
174+
else
175+
colormapfun = @winter;
176+
end
177+
178+
figure;
179+
subplot(2,2,1);
180+
imagesc(Sig_xx); colorbar; axis equal image;
181+
title('Correlation');
182+
subplot(2,2,2);
183+
imagesc(results.Sigma); colorbar; axis equal image;
184+
title('Conditional Corr.');
185+
subplot(2,2,3);
186+
imagesc(results.nCov); colorbar; axis equal image;
187+
title('Nuisance Factor');
188+
subplot(2,2,4);
189+
imagesc(options.corrfun(Xgsr)); colorbar; axis equal image;
190+
title('Correlation (NSR)');
191+
fullfile(options.outputdir,options.filename)
192+
print('-dpdf','-r150',fullfile(options.outputdir,options.filename)); pause(5)
193+
close all;
194+
end
195+
196+
197+
198+
end

‎src/+covariance/var_corr.m

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
function [R results] = var_corr(Sigma,options)
2+
%VAR_CORR - returns variance correlation decomposition of any symmetric, positive semi-definite matrix (usually covariance matrix)
3+
% Sigma = Gamma^{1/2} R Gamma^{1/2}
4+
% USAGE
5+
%
6+
% INPUTS
7+
%
8+
% OUTPUTS
9+
10+
if(~exist('options','var'))
11+
options = [];
12+
end
13+
14+
if(isfield(options,'verbose'))
15+
verbose=options.verbose;
16+
else
17+
options.verbose = true;
18+
verbose = options.verbose;
19+
end
20+
21+
import covariance.check_symposdef
22+
[issym isposdef results.symposdef] = check_symposdef(Sigma);
23+
try
24+
assert(issym>=1,'Not symmetric')
25+
catch me
26+
if(verbose)
27+
disp(me)
28+
results.symposdef.sym_err
29+
end
30+
end
31+
try
32+
assert(isposdef>=1,'Negative Eigenvalues')
33+
catch me
34+
if(verbose)
35+
disp(me)
36+
min(results.symposdef.eigs)
37+
end
38+
end
39+
40+
results.input.Sigma = Sigma;
41+
results.input.options = options;
42+
results.issym = issym;
43+
results.isposdef = isposdef;
44+
45+
if(issym==1)
46+
if(verbose)
47+
disp('Frob. Error: A-At')
48+
results.symposdef.sym_err
49+
disp('Symmetricizing input');
50+
end
51+
Sigma = (Sigma + Sigma')/2;
52+
end
53+
54+
p = length(Sigma);
55+
56+
gam = sqrt(diag(Sigma));
57+
d = 1./gam;
58+
59+
R = diag(d)*Sigma*diag(d);
60+
61+
results.var = d;
62+
results.corr = R;
63+
64+
end

‎src/+standardize/LICENSE.md

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
Copyright (C) 2017 Manjari Narayan
2+
3+
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
4+
5+
1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
6+
7+
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8+
9+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

0 commit comments

Comments
 (0)