Skip to content

Commit c5b8bd5

Browse files
XinhuiLiXinhuiLi
committed
Add MSIVA scripts for ISBI23 paper
1 parent 00415bb commit c5b8bd5

File tree

141 files changed

+17159
-2
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

141 files changed

+17159
-2
lines changed

‎.gitignore‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
# Other
2+
.vscode/
3+
run/
4+
model/__pycache__/
5+
model/__pycache__
6+
# Byte-compiled / optimized / DLL files
7+
__pycache__/
8+
*.py[cod]
9+
*$py.class
10+
11+
# C extensions
12+
*.so
13+
14+
# Distribution / packaging
15+
.Python
16+
build/
17+
develop-eggs/
18+
dist/
19+
downloads/
20+
eggs/
21+
.eggs/
22+
lib/
23+
lib64/
24+
parts/
25+
sdist/
26+
var/
27+
wheels/
28+
pip-wheel-metadata/
29+
share/python-wheels/
30+
*.egg-info/
31+
.installed.cfg
32+
*.egg
33+
MANIFEST
34+
35+
# PyInstaller
36+
# Usually these files are written by a python script from a template
37+
# before PyInstaller builds the exe, so as to inject date/other infos into it.
38+
*.manifest
39+
*.spec
40+
41+
# Installer logs
42+
pip-log.txt
43+
pip-delete-this-directory.txt
44+
45+
# Unit test / coverage reports
46+
htmlcov/
47+
.tox/
48+
.nox/
49+
.coverage
50+
.coverage.*
51+
.cache
52+
nosetests.xml
53+
coverage.xml
54+
*.cover
55+
*.py,cover
56+
.hypothesis/
57+
.pytest_cache/
58+
59+
# Translations
60+
*.mo
61+
*.pot
62+
63+
# Django stuff:
64+
*.log
65+
local_settings.py
66+
db.sqlite3
67+
db.sqlite3-journal
68+
69+
# Flask stuff:
70+
instance/
71+
.webassets-cache
72+
73+
# Scrapy stuff:
74+
.scrapy
75+
76+
# Sphinx documentation
77+
docs/_build/
78+
79+
# PyBuilder
80+
target/
81+
82+
# Jupyter Notebook
83+
.ipynb_checkpoints
84+
85+
# IPython
86+
profile_default/
87+
ipython_config.py
88+
89+
# pyenv
90+
.python-version
91+
92+
# pipenv
93+
# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
94+
# However, in case of collaboration, if having platform-specific dependencies or dependencies
95+
# having no cross-platform support, pipenv may install dependencies that don't work, or not
96+
# install all needed dependencies.
97+
#Pipfile.lock
98+
99+
# PEP 582; used by e.g. github.com/David-OConnor/pyflow
100+
__pypackages__/
101+
102+
# Celery stuff
103+
celerybeat-schedule
104+
celerybeat.pid
105+
106+
# SageMath parsed files
107+
*.sage.py
108+
109+
# Environments
110+
.env
111+
.venv
112+
env/
113+
venv/
114+
ENV/
115+
env.bak/
116+
venv.bak/
117+
118+
# Spyder project settings
119+
.spyderproject
120+
.spyproject
121+
122+
# Rope project settings
123+
.ropeproject
124+
125+
# mkdocs documentation
126+
/site
127+
128+
# mypy
129+
.mypy_cache/
130+
.dmypy.json
131+
dmypy.json
132+
133+
# Pyre type checker
134+
.pyre/
135+
.vscode/launch.json
136+
137+
# figure files
138+
*.pdf
139+
*.png
140+
*.p
141+
*.pt
142+
*.npy
143+
*.zip
144+
145+
.DS_Store

‎@MISAK/IVAfy.m‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
function Sold = IVAfy(O)
2+
3+
Sold = O.S;
4+
5+
S_ = cell(size(O.S));
6+
for mm = O.M
7+
ixs = [];
8+
for kk = find(O.nes)'
9+
ixs = [ixs find(O.S{mm}(kk,:) == 1, 1)];
10+
end
11+
S_{mm} = O.S{mm}(:,ixs);
12+
end
13+
O.updateCS(S_); % Update S and C
14+
15+
end

‎@MISAK/MISAK.m‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
classdef MISAK < handle
2+
% MISA using Kotz distribution and intrinsic source scaling
3+
properties
4+
C % Vector, number of components per dataset
5+
K % Number, number of non-empty subspaces
6+
M % Vector, indices of datasets to be used in the analysis
7+
N % Number, number of observations/measurements/samples
8+
S % Matrix, subspace assignment matrix (K x Cm)
9+
V % Vector, dimensionality of each dataset (number of sensors)
10+
W % Matrix, unmixing matrices (1 per dataset)
11+
X % Matrix, all datasets
12+
Y % Matrix, source estimates (network output)
13+
beta % Vector, Kotz beta parameter for each subspace
14+
d % Vector, dimensionality of each subspace
15+
d_k % Vector, dimensionality of each subspace per dataset
16+
eta % Vector, Kotz eta parameter for each subspace
17+
gradtype % String, type of gradient used: regular, relative, or natural
18+
sc % Logical, scale-control option: true or false
19+
lambda % Vector, Kotz lambda parameter for each subspace
20+
nes % Vector, non-empty subspace indexes
21+
nu % Vector, Kotz nu parameter for each subspace
22+
a % Vector, scaling factor to convers inverse cov. into inverse dispersion mtx.
23+
preX % Logical, preprocessing option: true, false, or empty
24+
end
25+
properties (Access = protected)
26+
ut % Container for utility static functions
27+
end
28+
methods
29+
function obj = MISAK(w0, M, S, X, beta, eta, lambda, gradtype, sc, preX)
30+
31+
%%% Test nu for all subspaces!!
32+
Sfull = full(sum([S{M}],2));
33+
Smask = Sfull~=0;
34+
if any(beta(Smask) <= 0), error('All beta parameters should be positive.'); end
35+
if ~isempty(lambda)
36+
if any(lambda(Smask) <= 0), error('All lambda parameters should be positive.'); end
37+
end
38+
if any(eta(Smask) <= ((2-Sfull(Smask))./2)), error('All eta parameters should be lagerer than (2-d)/2.'); end
39+
nu = (2*eta + Sfull - 2)./(2*beta);
40+
if any(nu(Smask) <= 0), error('All nu parameter derived from eta and d should be positive.'); end
41+
42+
obj.M = M; % Indices of datasets to be used in the analysis
43+
obj.S = S; % Subspace assignment matrix (K x Cm)
44+
obj.X = X; % All datasets
45+
obj.beta = beta; % Kotz beta parameter
46+
obj.eta = eta; % Kotz eta parameter
47+
obj.lambda = lambda; % Kotz lambda parameter
48+
obj.nu = nu; % Kotz nu parameter
49+
50+
if ~isempty(preX) && preX == true
51+
% Remove mean from all datasets:
52+
obj.preX.mean = cellfun(@(x) mean(x,2), obj.X, 'Un', 0);
53+
obj.X = cellfun(@(x,mx) bsxfun(@minus,x,mx), obj.X, obj.preX.mean, 'Un', 0);
54+
55+
% Standardize all datasets:
56+
obj.preX.std = cellfun(@(x) std(x,[],2), obj.X, 'Un', 0);
57+
obj.X = cellfun(@(sx,x) bsxfun(@times,1./sx,x), obj.preX.std, obj.X, 'Un', 0);
58+
59+
else obj.preX = false;
60+
end
61+
62+
obj.C = cellfun(@(s) size(s,2), obj.S); % Number of components per dataset
63+
obj.V = cellfun(@(x) size(x,1), obj.X); % Dimensionality of each dataset (number of sensors)
64+
65+
obj.ut = utils;
66+
67+
% Unmixing matrices (1 per dataset):
68+
obj.W = cell(1,max(obj.M));
69+
obj.W(obj.M) = obj.unstackW(w0);%, obj.M, obj.C(obj.M), obj.V(obj.M));
70+
71+
% Compute source estimates (network output):
72+
obj.Y = cell(1,max(obj.M));
73+
obj.Y(obj.M) = cellfun(@mtimes, obj.W(obj.M), obj.X(obj.M), 'Un', 0);
74+
75+
obj.N = size(obj.X{obj.M(1)},2);% Number of samples
76+
77+
obj.d = full(sum([obj.S{obj.M}],2)); % Dimensionality of each subspace
78+
obj.nes = obj.d~=0; % Non-empty subspace indexes
79+
obj.d_k = cellfun(@(s) full(sum(s,2)), obj.S,'Un',0);
80+
obj.auto_tune('lambda', lambda);
81+
82+
obj.a = (obj.lambda.^(-1./(obj.beta)) .* gamma(obj.nu + 1./obj.beta)) ./ ...
83+
(obj.d .* gamma(obj.nu));
84+
85+
obj.K = sum(obj.nes); % Number of non-empty subspaces
86+
if ~isempty(sc)
87+
updatesc(obj,sc); % Determines if scale-control is used: true or false
88+
else
89+
updatesc(obj,true); % Determines if scale-control is used: true or false
90+
end
91+
updategradtype(obj,gradtype); % Determines the type of gradient used: regular, relative, natural
92+
93+
end
94+
[J, gJ] = objective_(O) % Compute objective function value at current W
95+
[J, gJ] = objective_sc_(O) % Compute scale-control objective function value at current W
96+
[J, gJ] = objective(O,w) % Compute objective function value at provided W
97+
W = unstackW(O,w)
98+
w = stackW(O,W)
99+
update(O,S,M,b,l,e) % Subset the data
100+
Sold = IVAfy(O) % Turn problem into Unidimensional-type
101+
updateCS(O,S)
102+
updategradtype(O,gradtype) % Set gradient type
103+
combinatorial_optim(O,myM) % Solve combinatorial optimization
104+
w0 = greedysearch(O,myM)
105+
[w0, shuff] = sub_perm_analysis(O, w0)
106+
[shuff] = greedy_sub_perm_analysis(O)
107+
mISI = MISI(O,A) % Compute joint ISI of current W. A is provided.
108+
mmd = MMD(O,A)
109+
mmse = MMSE(O,Y)
110+
end
111+
methods (Access = private)
112+
% output = myFunc(obj,arg1,arg2)
113+
end
114+
methods (Static)
115+
% w = stackW(W)
116+
% W = unstackW(w,M,C,V)
117+
% Vt = myorth(Y)
118+
end
119+
end % End of classdef
120+
121+
% function myUtilityFcn
122+
% end

‎@MISAK/MISI.m‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
function [mISI] = MISI(O,A)
2+
3+
S = cellfun(@(s) full(s), O.S, 'Un', 0);
4+
mISI = O.ut.MISI(O.W,A,S);
5+
6+
end

‎@MISAK/MMD.m‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
function [mMD] = MMD(O,A)
2+
3+
mMD = O.ut.MMD(O.W,A,O.S);
4+
5+
end

‎@MISAK/MMSE.m‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
function [mMSE] = MMSE(O,Y)
2+
3+
mMSE = O.ut.MMSE(Y,O.Y,O.S);
4+
5+
end

‎@MISAK/auto_tune.m‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
function auto_tune(O, param, value)
2+
3+
if strcmpi(param,'eta')
4+
eta = value;
5+
if isempty(eta)
6+
eta = ones(length(O.nes),1);%ones(size(O.S{O.M(1)},1),1);
7+
end
8+
if length(eta) ~= length(O.nes)%size(O.S{O.M(1)},1)
9+
eta = eta(1).*ones(length(O.nes),1);%ones(size(O.S{O.M(1)},1),1);
10+
end
11+
min_val = (2-sum([O.S{O.M}],2))./2;
12+
if any(eta(O.nes) <= min_val(O.nes)), error('All eta parameters should be lagerer than (2-d)/2.'); end
13+
O.eta = eta; % Kotz eta parameter
14+
15+
return
16+
end
17+
18+
if strcmpi(param,'beta')
19+
beta = value;
20+
if isempty(beta)
21+
% beta = gammaln(pi)*3/4; % 0.617263173999417;
22+
% Get kurtosis-secific beta
23+
% 3 * (Formula 26) in this paper (3* is to scale to typical excess-kurtosis scale):
24+
% K. Zografos (2008), On Mardia�s and Song�s measures of kurtosis in elliptical distributions
25+
beta = zeros(length(O.nes),1);%zeros(size(O.S{O.M(1)},1),1);
26+
desired_kurtosis = 1.2; %1.591131166811146; % value that matches excess-kurtosis of logistic: ex-k = 1.2
27+
options = optimset('TolX',sqrt(eps));
28+
for kk = find(O.nes')
29+
%f = @(bt) 3*((((O.d(kk)^2)./(O.d(kk)*(O.d(kk)+2))) * ((gamma((2*O.eta(kk)+O.d(kk)+2)/(2*bt))*gamma((2*O.eta(kk)+O.d(kk)-2)/(2*bt)))/(gamma((2*O.eta(kk)+O.d(kk))/(2*bt))^2))) - 1) - desired_kurtosis;
30+
t1 = 2*log(O.d(kk))-log(O.d(kk))-log(O.d(kk)+2);
31+
g = @(bt) 3*(exp(t1 + (gammaln((2*O.eta(kk)+O.d(kk)+2)/(2*bt)) + gammaln((2*O.eta(kk)+O.d(kk)-2)/(2*bt)) - 2*gammaln((2*O.eta(kk)+O.d(kk))/(2*bt)))) - 1) - desired_kurtosis;
32+
beta(kk) = fzero(g,[7e4*sqrt(eps) 3], options);
33+
end
34+
end
35+
if length(beta) ~= length(O.nes)%size(O.S{O.M(1)},1)
36+
beta = beta(1).*ones(length(O.nes),1);%ones(size(O.S{O.M(1)},1),1);
37+
end
38+
if any(beta(O.nes) <= 0), error('All beta parameters should be positive.'); end
39+
O.beta = beta; % Kotz beta parameter
40+
41+
return
42+
end
43+
44+
if strcmpi(param,'lambda')
45+
lambda = value;
46+
if isempty(lambda)
47+
% lambda = pi/6; % 0.523541354951965;
48+
% Get variance-secific lambda
49+
% Using the definition of alpha in the MISA paper
50+
lambda = zeros(length(O.nes),1);%zeros(size(O.S{O.M(1)},1),1);
51+
desired_stddev = pi/sqrt(3); % value that matches stddev of logistic: std = pi/sqrt(3)
52+
% options = optimset('TolX',1e-100);
53+
for kk = find(O.nes')
54+
%f = @(lmb) ((lmb.^(-1./O.beta(kk)) .* gamma(O.nu(kk) + 1./O.beta(kk))) ./ (O.d(kk) .* gamma(O.nu(kk)))) - desired_stddev.^2;
55+
%g = @(lmb) gammaln(O.nu(kk) + 1./O.beta(kk)) - 1./O.beta(kk).*log(lmb) - log(O.d(kk)) - gammaln(O.nu(kk)) - 2*log(desired_stddev);
56+
%lambda(kk) = fzero(g,[0.1 100], options);
57+
lambda(kk) = exp(O.beta(kk) * (gammaln(O.nu(kk) + 1./O.beta(kk)) - log(O.d(kk)) - gammaln(O.nu(kk)) - 2*log(desired_stddev)));
58+
end
59+
end
60+
if length(lambda) ~= length(O.nes)%size(O.S{O.M(1)},1)
61+
lambda = lambda(1).*ones(length(O.nes),1);%ones(size(O.S{O.M(1)},1),1);
62+
end
63+
if any(lambda(O.nes) <= 0), error('All lambda parameters should be positive.'); end
64+
O.lambda = lambda; % Kotz lambda parameter
65+
66+
return
67+
end
68+
69+
end

0 commit comments

Comments
 (0)