Skip to content

Commit 75ef4e9

Browse files
FabienBouxFabienBoux
committed
Updated
1 parent 85df959 commit 75ef4e9

12 files changed

+218
-33
lines changed

‎Experiments/AliasingNoise.m

+4-3
Original file line numberDiff line numberDiff line change
@@ -219,9 +219,10 @@
219219

220220
disp(['Subsampling factor = ' num2str(nb_interleaves)])
221221

222-
SNR = 1./(3*nb_interleaves-2)*100;
223-
Irec = reshape(AddAliasingNoise(reshape(Iacq,[],size(Iacq,3)), SNR), size(Iacq,1),size(Iacq,2),size(Iacq,3));
224-
222+
SNR = 1./(3*nb_interleaves-2)*100;
223+
[Irec,real_snr] = AddAliasingNoise(reshape(Iacq,[],size(Iacq,3)), SNR);
224+
Irec = reshape(Irec, size(Iacq,1),size(Iacq,2),size(Iacq,3));
225+
225226
if S >= 430
226227
imgrec{nb_interleaves,1} = Irec(:,:,5);
227228
imgrec{nb_interleaves,2} = Irec(:,:,140);

‎Experiments/BoundaryBehavior_vMRF.m

+8-16
Original file line numberDiff line numberDiff line change
@@ -45,10 +45,7 @@
4545
% Quasi-random
4646
load(dicos{1})
4747
Xqmc = Dico.MRSignals(:,end/2+1:end) ./ Dico.MRSignals(:,1:end/2);
48-
% parfor i = 1:size(Xqmc,1)
49-
% Tmp(i,:) = interp1(Dico.Tacq(1:size(Xqmc,2)), Xqmc(i,:), echtime);
50-
% end
51-
% Xqmc = Tmp; Tmp = [];
48+
Xqmc = interp1(Dico.Tacq(1:size(Xqmc,2))', Xqmc', echtime)';
5249
Yqmc = Dico.Parameters.Par;
5350

5451
Xqmc(Dico.Parameters.Par(:,4)== 1e-6,:) = [];
@@ -69,10 +66,7 @@
6966
% Grid
7067
load(dicos{2})
7168
Xgrid = Dico.MRSignals(:,end/2+1:end) ./ Dico.MRSignals(:,1:end/2);
72-
% parfor i = 1:size(Xgrid,1)
73-
% Tmp(i,:) = interp1(Dico.Tacq(1:size(Xgrid,2)), Xgrid(i,:), echtime);
74-
% end
75-
% Xgrid = Tmp; Tmp = [];
69+
Xgrid = interp1(Dico.Tacq(1:size(Xgrid,2))', Xgrid', echtime)';
7670
Ygrid = Dico.Parameters.Par;
7771

7872
Xgrid(Dico.Parameters.Par(:,4)== 1e-6,:) = [];
@@ -93,10 +87,7 @@
9387
% Random - test signals
9488
load(dicos{3})
9589
Xtest_ = Dico.MRSignals(:,end/2+1:end) ./ Dico.MRSignals(:,1:end/2);
96-
% parfor i = 1:size(Xtest_,1)
97-
% Tmp(i,:) = interp1(Dico.Tacq(1:size(Xtest_,2)), Xtest_(i,:), echtime);
98-
% end
99-
% Xtest_ = Tmp; Tmp = [];
90+
Xtest_ = interp1(Dico.Tacq(1:size(Xtest_,2))', Xtest_', echtime)';
10091
Ytest_ = Dico.Parameters.Par;
10192

10293
Xtest_(Dico.Parameters.Par(:,4)== 1e-6,:) = [];
@@ -115,8 +106,8 @@
115106

116107
%% Reduce datasets to keep only signals with StO2 around 0.7
117108

118-
bds_sto2 = [0.65 0.75];
119-
bds_sto2 = [0.6 0.8];
109+
% bds_sto2 = [0.65 0.75];
110+
bds_sto2 = [0 1];
120111

121112
v = (Ygrid(:,3) <= bds_sto2(1) | Ygrid(:,3) >= bds_sto2(2));
122113
Xgrid(v,:) = [];
@@ -148,8 +139,9 @@
148139
13e-2 23e-6 ;
149140
16e-2 24e-6 ;
150141
19e-2 25e-6 ;
151-
22e-2 26e-6 ;
152-
... %then, boundary entries
142+
22e-2 26e-6 ;];
143+
ext = [ext; ext + ext.*0.1.*rand(size(ext)) + ext.*0.2.*rand(size(ext))];
144+
ext = [ ext;... %then, boundary entries
153145
max(Yqmc(:,1)) 5e-6;
154146
max(Yqmc(:,1)) max(Yqmc(:,2));
155147
15e-2 max(Yqmc(:,2));

‎Experiments/Kimpact_vMRF.m

+157
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,157 @@
1+
2+
%% Description
3+
%
4+
% TODO
5+
%
6+
% Fabien Boux - 11/2020
7+
8+
Init
9+
disp(['Running experiment ' mfilename '.m'])
10+
11+
12+
%% Setting
13+
14+
% Execution settings
15+
verbose = 1;
16+
backup = 1;
17+
18+
% Experiment settings
19+
nb_test_signals = 5e4;
20+
nb_signals = [2e3 1e4];
21+
22+
K = 3:3:200;
23+
constraints = {'d*', 'i*'};
24+
25+
% Regression settings
26+
Model.Lw = 0;
27+
Parameters.cstr.Gammat = '';
28+
Parameters.cstr.Gammaw = '';
29+
Model.maxiter = 250;
30+
Model.snrtrain = 60;
31+
32+
33+
%% Load data set
34+
35+
%dictionary filenames
36+
dicos = {'inputs/3DvMRF/qMC.mat'}; %qMC
37+
38+
% Load echotimes
39+
echtime = load('inputs/echoetime.mat');
40+
echtime = echtime.t;
41+
42+
% Quasi-random
43+
load(dicos{1})
44+
Xqmc = Dico.MRSignals(:,end/2+1:end) ./ Dico.MRSignals(:,1:end/2);
45+
Xqmc = interp1(Dico.Tacq(1:size(Xqmc,2))', Xqmc', echtime)';
46+
Yqmc = Dico.Parameters.Par;
47+
48+
Xqmc(Dico.Parameters.Par(:,4)== 1e-6,:) = [];
49+
Yqmc(Dico.Parameters.Par(:,4)== 1e-6,:) = [];
50+
51+
Yqmc = Yqmc(:, [2 3 1]);
52+
53+
[rows, ~] = find(~isnan(Xqmc));
54+
r = unique([rows; find(any(isnan(Yqmc),2) == 0)]);
55+
Xqmc = Xqmc(r,:);
56+
Yqmc = Yqmc(r,:);
57+
58+
Xtest = Xqmc(nb_signals(end)+1:end,:);
59+
Ytest = Yqmc(nb_signals(end)+1:end,:);
60+
61+
% Add noise
62+
Xqmc = AddNoise(Xqmc, Model.snrtrain);
63+
Xtest = AddNoise(Xtest, 50);
64+
65+
66+
%% Creating data
67+
68+
% Err = struct(length(nb_params),length(constraints));
69+
% Bic = struct(length(nb_params),length(constraints));
70+
71+
for c1 = 1:length(constraints)
72+
for c2 = 1:length(nb_signals)
73+
74+
% Choose number of parameters and constrint on covariance matrices
75+
Model.cstr.Sigma = constraints{c1};
76+
77+
% Generate training dataset
78+
Xtrain = Xqmc(1:nb_signals(c2),:);
79+
Ytrain = Yqmc(1:nb_signals(c2),:);
80+
81+
% Init
82+
Rmse = nan(length(K), size(Ytrain,2));
83+
bic = nan(1,length(K));
84+
[N,Lt] = size(Ytrain);
85+
Lw = Model.Lw;
86+
D = size(Xtrain,2);
87+
88+
% Start evaluating K impact
89+
parfor k = 1:length(K)
90+
disp(k)
91+
92+
try
93+
v = randi(size(Xtest,1),nb_test_signals,1);
94+
95+
[theta,~,ll] = EstimateInverseFunction(Ytrain, Xtrain, K(k), Model.Lw, Model.maxiter, Model.cstr, 0);
96+
Ypredict = EstimateParametersFromModel(Xtest(v,:), theta, 0);
97+
98+
Rmse(k,:) = EvaluateEstimation(Ytest(v,:), Ypredict);
99+
100+
switch Model.cstr.Sigma
101+
case 'i*'
102+
M = K(k)* (D*(Lw+Lt+1) + Lt*(Lt+3)/2 + 1);
103+
case 'i'
104+
M = K(k)* (D*(Lw+Lt+1) + Lt*(Lt+3)/2 + 1) + K(k)-1;
105+
case 'd*'
106+
M = K(k)* (D*(Lw+Lt+2) + Lt*(Lt+3)/2) + D + K(k)-1;
107+
case 'd'
108+
M = K(k)* (D*(Lw+Lt+2) + Lt*(Lt+3)/2) + K(k)-1;
109+
otherwise
110+
M = K(k)* (D*(Lw+Lt+D+1) + Lt*(Lt+3)/2) + K(k)-1;
111+
end
112+
113+
bic(k) = -2*ll + M*log(N);
114+
115+
catch
116+
warning(['Iteration ' Model.cstr.Sigma '/' num2str(K(k)) ' failed'])
117+
end
118+
end
119+
120+
Err{c2,c1} = Rmse;
121+
Bic{c2,c1} = bic;
122+
end
123+
end
124+
125+
126+
%% Backup
127+
128+
if backup == 1
129+
clear tmp* Dico* Estim X* Y*
130+
save(['temp/' 'Kimpact_vMRF'])
131+
end
132+
133+
134+
%% Disp
135+
136+
fig = figure;
137+
hold on
138+
for p = 1:3
139+
subplot(1,3,p)
140+
hold on
141+
142+
for c2 = 1:length(nb_signals)
143+
plot(K, mean(Err{c2,1}(:,p),2), '.-', 'linewidth',1.5);
144+
%plot(K, mean(Dic{c2,c1},2), '.-', 'linewidth',1.5);
145+
end
146+
end
147+
ylabel('RMSE')
148+
%ylabel('BIC')
149+
xlabel('K')
150+
% legend([repelem('P = ',length(nb_params),1) num2str(nb_params')])
151+
152+
153+
%%
154+
155+
if backup == 1
156+
savefig(fig, ['figures/' 'Kimpact_vMRF'])
157+
end

‎Experiments/ParameterSpaceSampling.m

+8-4
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,7 @@
3737
snr_test = inf;
3838

3939
% Method settings
40-
methods = {'DBM', 'DB-DL'}; %{'DBM', 'DB-SL'};
40+
methods = {'DBM', 'DB-DL', 'DB-SL'};
4141
sampling_strategies = {'Grid', 'Random', 'qRandom'};
4242
Model.snrtrain = inf;
4343
Model.Exec = 'cpu';
@@ -160,7 +160,7 @@
160160
if exist(filename,'file')
161161
load(filename,'mRMSE_DBM','mRMSE_DBSL','mRMSE_DBDL','nb_param','nb_train_signals');
162162
DBM_rmse{i} = mRMSE_DBM;
163-
DBL_rmse{i} = mRMSE_DBDL;
163+
DBL_rmse{i} = mRMSE_DBSL;
164164
else
165165
DBM_rmse{i} = [];
166166
DBL_rmse{i} = [];
@@ -179,7 +179,7 @@
179179
for i = [1 4 7]
180180
c = c + 1;
181181

182-
dat = DBL_rmse{i}';
182+
dat = DBM_rmse{i}';
183183
h(c) = subplot(1,3,c);
184184

185185
boxplot(dat,'symbol','')
@@ -193,9 +193,13 @@
193193

194194
title([titles{c} ' ' num2str(title_nb_param(i)) ' parameters'])
195195

196-
% m = mean(dat);
196+
m = mean(dat);
197+
s = std(dat);
197198
% disp(['qRand / Grid: ' num2str(1 - m(3) / m(1))])
198199
% disp(['qRand / Rand: ' num2str(1 - m(3) / m(2))])
200+
disp(['P = ' num2str(title_nb_param(i))])
201+
disp(1e3*m)
202+
disp(1e3*s)
199203
end
200204

201205
linkaxes(h,'y')

‎Experiments/vMRFperformance.m

+39-8
Original file line numberDiff line numberDiff line change
@@ -19,12 +19,12 @@
1919
nb_train_signals = 150500;
2020
nb_test_signals = 100000;
2121

22-
snr_test = 100;
22+
snr_test = 50;
2323
snr_train = 60;
2424

25-
dicos = {{'inputs/new/qMC.mat', ' '},... %qMC
26-
{'inputs/new/Regular.mat', ' '},... %grid
27-
{'inputs/new/qMC.mat', ' '},... %test signals
25+
dicos = {{'inputs/3DvMRF/qMC.mat', ' '},... %qMC
26+
{'inputs/3DvMRF/Regular.mat', ' '},... %grid
27+
{'inputs/3DvMRF/qMC.mat', ' '},... %test signals
2828
};
2929
use_recomputed_parameter = 0;
3030

@@ -61,8 +61,9 @@
6161
% Quasi-random
6262
load(dicos{1}{s})
6363
if use_recomputed_parameter == 1
64-
Dico.Parameters.Par(:,1) = Dico.Parameters.Par(:,7);
6564
Dico.Parameters.Par(:,2) = Dico.Parameters.Par(:,8);
65+
Dico.Parameters.Par(:,3) = Dico.Parameters.Par(:,9);
66+
Dico.Parameters.Par(:,4) = Dico.Parameters.Par(:,10);
6667
end
6768
Xqmc = Dico.MRSignals(:,end/2+1:end) ./ Dico.MRSignals(:,1:end/2);
6869
parfor i = 1:size(Xqmc,1)
@@ -85,8 +86,9 @@
8586
% Regular
8687
load(dicos{2}{s})
8788
if use_recomputed_parameter == 1
88-
Dico.Parameters.Par(:,1) = Dico.Parameters.Par(:,7);
8989
Dico.Parameters.Par(:,2) = Dico.Parameters.Par(:,8);
90+
Dico.Parameters.Par(:,3) = Dico.Parameters.Par(:,9);
91+
Dico.Parameters.Par(:,4) = Dico.Parameters.Par(:,10);
9092
end
9193
Xgrid = Dico.MRSignals(:,end/2+1:end) ./ Dico.MRSignals(:,1:end/2);
9294
parfor i = 1:size(Xgrid,1)
@@ -109,8 +111,9 @@
109111
% Picking test signals
110112
load(dicos{3}{s})
111113
if use_recomputed_parameter == 1
112-
Dico.Parameters.Par(:,1) = Dico.Parameters.Par(:,7);
113114
Dico.Parameters.Par(:,2) = Dico.Parameters.Par(:,8);
115+
Dico.Parameters.Par(:,3) = Dico.Parameters.Par(:,9);
116+
Dico.Parameters.Par(:,4) = Dico.Parameters.Par(:,10);
114117
end
115118
Xtest_ = Dico.MRSignals;
116119
parfor i = 1:size(Xtest_,1)
@@ -131,7 +134,7 @@
131134

132135
%% DBL method training
133136

134-
Model_.K = 50;
137+
Model_.K = 200;
135138
Model_.Lw = 0;
136139

137140
r = randperm(size(Xqmc,1));
@@ -363,3 +366,31 @@
363366
if backup == 1
364367
savefig(fig, ['figures/' 'vMRFperformance'])
365368
end
369+
370+
371+
%% Values
372+
373+
disp('%%% BVf RMSE %%')
374+
disp(['CEF: ' num2str(nanmean(1e2*reshape(Rmse_anlt(:,:,:,1), 1,[])), '%.2f') ' %'])
375+
disp(['DBM: ' num2str(nanmean(1e2*reshape(Rmse_grid(:,:,:,1), 1,[])), '%.2f') ' %'])
376+
disp(['DBDL: ' num2str(nanmean(1e2*reshape(Rmse_nn(:,:,:,1), 1,[])), '%.2f') ' %'])
377+
disp(['DBSL: ' num2str(nanmean(1e2*reshape(Rmse_gllim(:,:,:,1), 1,[])), '%.2f') ' %'])
378+
disp(['DBSL/CEF: ' num2str(nanmean(1e2*reshape(Rmse_anlt(:,:,:,1) - Rmse_gllim(:,:,:,1), 1,[])), '%.2f') ' %'])
379+
disp(['DBSL/DBM: ' num2str(nanmean(1e2*reshape(Rmse_grid(:,:,:,1) - Rmse_gllim(:,:,:,1), 1,[])), '%.2f') ' %'])
380+
381+
disp(' ')
382+
disp('%% VSI RMSE %%')
383+
disp(['CEF: ' num2str(nanmean(1e6*reshape(Rmse_anlt(:,:,:,2), 1,[])), '%.2f') ' um'])
384+
disp(['DBM: ' num2str(nanmean(1e6*reshape(Rmse_grid(:,:,:,2), 1,[])), '%.2f') ' um'])
385+
disp(['DBDL: ' num2str(nanmean(1e6*reshape(Rmse_nn(:,:,:,2), 1,[])), '%.2f') ' um'])
386+
disp(['DBSL: ' num2str(nanmean(1e6*reshape(Rmse_gllim(:,:,:,2), 1,[])), '%.2f') ' um'])
387+
disp(['DBSL/CEF: ' num2str(nanmean(1e6*reshape(Rmse_anlt(:,:,:,2) - Rmse_gllim(:,:,:,2), 1,[])), '%.2f') ' um'])
388+
disp(['DBSL/DBM: ' num2str(nanmean(1e6*reshape(Rmse_grid(:,:,:,2) - Rmse_gllim(:,:,:,2), 1,[])), '%.2f') ' um'])
389+
390+
disp(' ')
391+
disp('%% CI %%')
392+
disp(['Max diff RMSE - CI: BVf= ' num2str(max(reshape(1e2*abs(Rmse_gllim(:,:,:,1) - CI(:,:,:,1)), 1,[])), '%.2f')...
393+
' %, VSI= ' num2str(max(reshape(1e6*abs(Rmse_gllim(:,:,:,2) - CI(:,:,:,2)), 1,[])), '%.2f') ' um'])
394+
disp(['Mean diff RMSE - CI: BVf= ' num2str(nanmean(reshape(1e2*abs(Rmse_gllim(:,:,:,1) - CI(:,:,:,1)), 1,[])), '%.2f')...
395+
' %, VSI= ' num2str(nanmean(reshape(1e6*abs(Rmse_gllim(:,:,:,2) - CI(:,:,:,2)), 1,[])), '%.2f') ' um'])
396+

‎figures/BoundaryBehaviour_vMRF_1.fig

2.02 MB
Binary file not shown.

‎figures/BoundaryBehaviour_vMRF_2.fig

5.88 KB
Binary file not shown.

‎figures/Kimpact_vMRF.fig

62.9 KB
Binary file not shown.

‎figures/toy_signals_FT.fig

101 KB
Binary file not shown.

‎figures/vMRFperformance.fig

4.19 KB
Binary file not shown.

‎tools/GLLiM-api/AddAliasingNoise.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616

1717
narginchk(2, 2);
1818

19-
if ~any(imag(X) ~= 0)
19+
if ~any(imag(X(:)) ~= 0)
2020
X_noisy = X + abs(X) .*(randn(size(X)) ./snr);
2121

2222
real_snr = max(X,[],2) ./ std(X - X_noisy, [],2);

‎tools/GLLiM-api/AddNoise.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
narginchk(2, 2);
1616

1717
% real signals
18-
if ~any(imag(X) ~= 0)
18+
if ~any(imag(X(:)) ~= 0)
1919
X_noisy = abs(X + randn(size(X)) .* repmat(max(abs(X),[],2) ./ snr, 1,size(X,2)));
2020

2121
real_snr = max(X,[],2) ./ std(X - X_noisy, [],2);

0 commit comments

Comments
 (0)