diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..14b1a71bff7ac236bbb42ee007315ab7c1eb4276
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 Max Planck Institute Group Anke Henning
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..0c39dcb500b456f03893cc6899fc1f20d547467f
--- /dev/null
+++ b/README.md
@@ -0,0 +1,54 @@
+## Description
+This library provides extensive functionality to reconstruct magnetic resonance spectroscopy data.
+These include:
+1. preprocess raw data including:
+ - metabolite cycling
+ - frequency alignment both in time and frequency domain
+ - Averaging
+ - Eddy current correction
+ - coil combination
+ - HSVD water suppression
+ - etc. See functions such as `reconstruct`, `reconstruct_all_...`
+2. call LCModel spectral fitting via ssh
+3. calculate mmol/kg concentrations
+4. all above functionalities are adapted to:
+ - upfield spectra
+ - downfield spectra
+ - macromolecular spectra
+ - macromolecular spectra fitted to amino acids models
+ - creating and evaluating simulated spectra
+ - functional MRS data
+
+## Dependencies
+
+For reading in raw scanner data files, the `mapVBVD` function by Philipp Ehses was used. This and other functions to read in raw scanner data can be found in the Gannet2.0 repository: https://github.com/cjohnevans/Gannet2.0
+Additional MATLAB packages used by the code include:
+https://de.mathworks.com/matlabcentral/fileexchange/27999-ssh-from-matlab-updated-sftp-scp
+https://de.mathworks.com/matlabcentral/fileexchange/19729-passwordentrydialog
+https://de.mathworks.com/matlabcentral/fileexchange/27418-fdr_bh
+https://de.mathworks.com/matlabcentral/fileexchange/28303-bonferroni-holm-correction-for-multiple-comparisons
+https://de.mathworks.com/matlabcentral/fileexchange/71052-blandaltmanplot
+https://de.mathworks.com/matlabcentral/fileexchange/13634-axescoord2figurecoord
+https://de.mathworks.com/matlabcentral/fileexchange/5782-eps2pdf
+
+## Publications
+
+Following noteable publications used this library:
+
+Giapitzakis, I. A., Shao, T., Avdievich, N., Mekle, R., Kreis, R., & Henning, A. (2018). Metaboliteâ€cycled STEAM and semiâ€LASER localization for MR spectroscopy of the human brain at 9.4 T. Magnetic resonance in medicine, 79(4), 1841-1850.
+
+Giapitzakis, I. A., Borbath, T., Muraliâ€Manohar, S., Avdievich, N., & Henning, A. (2019). Investigation of the influence of macromolecules and spline baseline in the fitting model of human brain spectra at 9.4 T. Magnetic resonance in medicine, 81(2), 746-758.
+
+Giapitzakis, I. A., Avdievich, N., & Henning, A. (2018). Characterization of macromolecular baseline of human brain using metabolite cycled semiâ€LASER at 9.4 T. Magnetic resonance in medicine, 80(2), 462-473.
+
+Avdievich, N. I., Giapitzakis, I. A., Pfrommer, A., Borbath, T., & Henning, A. (2018). Combination of surface and ‘vertical’loop elements improves receive performance of a human head transceiver array at 9.4 T. NMR in Biomedicine, 31(2), e3878.
+
+Muraliâ€Manohar, S., Borbath, T., Wright, A. M., Soher, B., Mekle, R., & Henning, A. (2020). T2 relaxation times of macromolecules and metabolites in the human brain at 9.4 T. Magnetic resonance in medicine, 84(2), 542-558.
+
+Muraliâ€Manohar, S., Wright, A. M., Borbath, T., Avdievich, N. I., & Henning, A. (2021). A novel method to measure T1â€relaxation times of macromolecules and quantification of the macromolecular resonances. Magnetic resonance in medicine, 85(2), 601-614.
+
+Borbath, T., Muraliâ€Manohar, S., Wright, A. M., & Henning, A. (2021). In vivo characterization of downfield peaks at 9.4 T: T2 relaxation times, quantification, pH estimation, and assignments. Magnetic resonance in medicine, 85(2), 587-600.
+
+Dorst, J., Ruhm, L., Avdievich, N., Bogner, W., & Henning, A. (2021). Comparison of four 31P singleâ€voxel MRS sequences in the human brain at 9.4 T. Magnetic resonance in medicine, 85(6), 3010-3026.
+
+Borbath, T., Muraliâ€Manohar, S., Dorst, J., Wright, A. M., & Henning, A. (2021). ProFitâ€1D—A 1D fitting software and openâ€source validation data sets. Magnetic resonance in medicine
diff --git a/algorithms/ExportLCModelBasis.m b/algorithms/ExportLCModelBasis.m
index 7ecfa191c6d52222692a21318e7529d971010c1a..6e1c91421dfca41c2bda559a2b17dcbc640cd6d0 100644
--- a/algorithms/ExportLCModelBasis.m
+++ b/algorithms/ExportLCModelBasis.m
@@ -51,7 +51,7 @@ flagY = false;
switch choice
case 'Yes'
flagY = true;
- time = linspace(0,dwellTimeMs*NSamples,NSamples);
+ time = linspace(0,dwellTimeMs*NSamples,NSamples); % calculation in milliseconds!
FWHM = 3; %Hz
T2 = 1/(pi*FWHM)*1e+3; %ms
% Mo = 4.9*(1+T2)/T2; %in order the peak to be 4.9 in the FD-real part
diff --git a/algorithms/ImportLCModelBasis.m b/algorithms/ImportLCModelBasis.m
index 4bd8cfd3522a65a09bce48e945d7f6af5f12450c..4efb06c8924b4daebc49cdbf190fc95e537d2212 100644
--- a/algorithms/ImportLCModelBasis.m
+++ b/algorithms/ImportLCModelBasis.m
@@ -12,9 +12,9 @@ if ~(strcmp(filepath(end),'\') || strcmp(filepath(end),'/'))
filepath = [filepath '/'];
end
-% if ~contains(filename,'.raw', 'IgnoreCase',true)
-% filename = strcat(filename, '.RAW');
-% end
+if ~contains(filename,'.raw', 'IgnoreCase',true)
+ filename = strcat(filename, '.RAW');
+end
%open the .RAW data
filename = strcat(filepath, filename);
%
diff --git a/createTestData/createSimulatedSpectra.m b/createTestData/createSimulatedSpectra.m
index a287c53efcb9c4dbab17c43db62c213c2652453a..171c4bbb1d8a05f348991832638174e232312bca 100644
--- a/createTestData/createSimulatedSpectra.m
+++ b/createTestData/createSimulatedSpectra.m
@@ -15,7 +15,7 @@ TE = 24 * 1e-3;
numMetabolites = length(metabolites);
saveData = false;
-saveFigures = true;
+saveFigures = false;
saveRandNumbers = false; %should generally be kept false to have the same random numbers consistently
%% path to the export of the spectra
@@ -35,19 +35,17 @@ pc1 = [-17.5:2.5:17.5]; %grad/ppm
% pc1 = [-30:5:30]; %grad/ppm
df2_global = [-17.5:2.5:17.5]; % Hz : Global and local shifts
df2_local_scale = [0.25 0.25 0.5 0.5 0.75 0.75 1 1 1.5 1.5 2 2 3 3 5];
-df2_local_std = 3;% 3 Hz coresponds to ~0.0075 ppm.
+df2_local_std = 3;% 3 Hz coresponds to ~0.0075 ppm.
% This is the maximal shift, which can occur by temperature between
% moieties (Wermter, MRMP 2017), assuming up to 10°C temperature change.
-% pH changes should have a similar order of magnitude.
+% pH changes should have a similar order of magnitude.
em_local_scale = [0.2 0.2 0.4 0.4 0.6 0.6 0.8 0.8 1 1 1.2 1.2 1.5 1.5 2];
-% when truncating at 200 ms like used in real spectra, the NAA SNR of the defaultSummedFid corresponds to
-% Chosen NoiseLevel -> SNR(NAA)
-% 1 -> SNR = 207; 3 -> SNR = 158; 7 -> SNR = 110; 15 -> SNR = 44; 25 -> SNR = 13
-noiseLevel = [1:25];
+%check the resulting SNR levels from noiseLevel in snr_NAA_noise in the noise simulation
+noiseLevel = [1:25];
numberDfPermutations = length(df2_global);
-if saveRandNumbers
+if saveRandNumbers
%generally you should not enter this if. It creates random numbers which can be used for calculating the concentrations
df2_local_RandNumbers = rand(numberDfPermutations,numMetabolites);
save([dataExportPathBase, 'df2_local_RandNumbers.mat'], 'df2_local_RandNumbers');
@@ -57,7 +55,7 @@ em_local_RandNumbers = df2_local_RandNumbers; %use the same random numbers for s
%%
metaboliteFileNames = metabolites;
-% replace the weirdos
+% replace the weirdos
metaboliteFileNames = strrep(metaboliteFileNames, 'tCr(CH2)', 'Cr_singlet_tot_3.925');
metaboliteFileNames = strrep(metaboliteFileNames, 'tCr(CH3)', 'Cr_singlet_tot_3.028');
metaboliteFileNames = strrep(metaboliteFileNames, 'NAA(CH3)', 'NAA_2');
@@ -89,7 +87,7 @@ for indexMetabolite = 1:numMetabolites
[fid, bandwidth, scanFrequency, metaboliteName, singletPresent] = ...
ImportLCModelBasis(filePathMetabolite, fileName, plotBasis, '1H');
if strcmp(metabolites{indexMetabolite},'NAA_ac')
-% fid = fid * 3;
+ % fid = fid * 3;
end
if strcmp(metabolites{indexMetabolite},'MMB')
fids{indexMetabolite}(1,:)= fid(1:4096); % stupid change to make all metabolites 4096 long!
@@ -113,6 +111,7 @@ phase_1_f2_mx = ppmVector-phasePivot_ppm;
ppmMask = (ppmVector>0.6) & (ppmVector<4.1);
+legendLabels = cell(3,1);
%% set default values for the parameters
default_pc0 = 0;
@@ -131,33 +130,37 @@ current_pc1 = default_pc1;
current_df2 = default_df2;
current_gm = default_gm;
current_em_std_factor = default_em_std_factor;
-conc_std_factor = zeros(1,numMetabolites);
+conc_std_factor = zeros(1,numMetabolites);
% generate the default spectrum
-[defaultSummedFid, defaultSummedFidWithNoise, defaultCurrent_em, defaultCurrentConcentrationRm] = ...
+[defaultSummedFid, defaultSummedFidWithNoise, defaultCurrent_em, defaultCurrentConcentrationRm, defaultStats] = ...
createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, 'meanConc', [], dwellTime * 1e3, scanFrequency, dataExportPathBase, water_ppm, scalingFactorRefPeak);
+plot(ppmVector,real(fftshift(fft(defaultSummedFid))));
+% figure(1);
+% subplot(5,4,8+7:8+8)
figure(1);
-subplot(5,4,8+7:8+8)
+subplot(3,4,9:10)%Revision 1
% figure(2);
% subplot(3,4,7:8)
-plot(ppmVector,real(fftshift(fft(defaultSummedFid))));
-
+
%% generate the series of different concentrations
stdLimConc = 3.5;
numberConcPermutations = 100;
concentrationsRm = zeros(numberConcPermutations, numMetabolites);
paramKeyword = 'conc';
-if saveRandNumbers
+titleStr = '{\boldmath$c$}';
+if saveRandNumbers
%generally you should not enter this if. It creates random numbers which can be used for calculating the concentrations
concentrationRandNumbers = rand(numberConcPermutations,numMetabolites);
save([dataExportPathBase, 'concentrationRandNumbers.mat'], 'concentrationRandNumbers');
end
load([dataExportPathBase, 'concentrationRandNumbers.mat'], 'concentrationRandNumbers');
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
+snr_NAA_conc = zeros(numberConcPermutations,1);
for indexParam = 1:numberConcPermutations
for indexMetabolite = 1:numMetabolites
% generate a random number between [-stdLim; stdLim]
@@ -176,133 +179,153 @@ for indexParam = 1:numberConcPermutations
end
conc_std_factor(indexMetabolite) = stdRand;
end
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
+ snr_NAA_conc(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
- plot(ppmVector,real(fftshift(fft(summedFid))));
+ plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
% ImportLCModelBasis(dataExportPath, fileName, true, '1H');
end
-titleStr = '{\boldmath$c$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-for indexMetabolite = 1:numMetabolites
- figure;
- scatter(concentrationsRm(:,indexMetabolite), concentrationsRm(:,indexMetabolite));
- title(metabolites(indexMetabolite))
-end
+% for indexMetabolite = 1:numMetabolites
+% figure;
+% scatter(concentrationsRm(:,indexMetabolite), concentrationsRm(:,indexMetabolite));
+% title(metabolites(indexMetabolite))
+% end
%% #############################################################################################
% generate the series of different parameters, same concentrations
% #############################################################################################
conc_std_factor = zeros(1,numMetabolites); % reset value
-%% gauss
+%% gauss
paramKeyword = 'gauss';
+titleStr = '$$\nu_{g}$$';
+units = '(Hz)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(gaussianValues), numMetabolites);
-figure(1);
-subplot(5,4,5:6)
+% figure(1);
+% subplot(5,4,5:6)
+figure(2);
+subplot(2,4,1:2)%Revision 1
% subplot(2,4,5:6)
+snr_NAA_Gauss = zeros(length(gaussianValues),1);
for indexParam = 1:length(gaussianValues)
current_gm = gaussianValues(indexParam);
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
-
+
+ snr_NAA_Gauss(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFid))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, false, legendLabels, indexParam, gaussianValues, '$\rm\nu_{g}$', units);
end
% reset parameter
current_gm = default_gm;
-
-titleStr = '$$\nu_{g}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
+legend(legendLabels, 'Interpreter', 'latex');
+plotArrowsPlotOffsets()
+%%
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-%% pc0 - zero order phase
+%% pc0 - zero order phase
paramKeyword = 'pc0';
+titleStr = '$$\varphi_{0}$$';
+units = '(degrees)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(pc0), numMetabolites);
figure(1);
-subplot(5,4,1:2)
+% subplot(5,4,1:2)
+subplot(3,4,1:2) %Revision 1
% subplot(2,4,1:2)
+snr_NAA_pc0 = zeros(length(pc0),1);
for indexParam = 1:length(pc0)
current_pc0 = pc0(indexParam);
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
-
+
+ snr_NAA_pc0(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFid))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, pc0, titleStr, units);
end
% reset parameter
current_pc0 = default_pc0;
-titleStr = '$$\varphi_{0}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-%% pc1 - first order phase
+%% pc1 - first order phase
paramKeyword = 'pc1';
+titleStr = '$$\varphi_{1}$$';
+units = '(degrees/ppm)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(pc1), numMetabolites);
figure(1);
-subplot(5,4,3:4)
+% subplot(5,4,3:4)
+subplot(3,4,3:4)%Revision 1
% subplot(2,4,3:4)
+snr_NAA_pc1 = zeros(length(pc1),1);
for indexParam = 1:length(pc1)
current_pc1 = pc1(indexParam);
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
-
+
+ snr_NAA_pc1(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFid))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, pc1, titleStr, units);
end
% reset parameter
current_pc1 = default_pc1;
-titleStr = '$$\varphi_{1}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-%% df2 local - frequency shifts
+%% df2 local - frequency shifts
paramKeyword = 'df2_local';
+titleStr = '{\boldmath$\omega_{local}$}';
+units = '(Hz)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(df2_local_scale), numMetabolites);
+% figure(1);
+% subplot(5,4,8+1:8+2)
figure(1);
-subplot(5,4,8+1:8+2)
+subplot(3,4,7:8)%Revision 1
% figure(2);
% subplot(3,4,1:2)
+snr_NAA_df2_local = zeros(length(df2_local_scale),1);
for indexParam = 1:length(df2_local_scale)
-
+
for indexMetabolite = 1:numMetabolites
% generate a random number between [-stdLim; stdLim]
stdRand = 2 * df2_local_std * df2_local_RandNumbers(indexParam,indexMetabolite) - df2_local_std;
@@ -316,67 +339,77 @@ for indexParam = 1:length(df2_local_scale)
case 'NAA(CH2)'
stdRand = NaaStdRand;
otherwise
- % nothing to do
+ % nothing to do
end
current_df2(indexMetabolite) = df2_local_scale(indexParam) * stdRand;
end
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, water_ppm, scalingFactorRefPeak);
-
+
+ snr_NAA_df2_local(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFid))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, df2_local_scale, titleStr, units);
end
% reset parameter
current_df2 = default_df2;
-titleStr = '{\boldmath$\omega_{local}$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-%% df2 global - frequency shifts
+%% df2 global - frequency shifts
paramKeyword = 'df2_global';
+titleStr = '$$\omega_{global}$$';
+units = ('Hz');
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(df2_global), numMetabolites);
+% figure(1);
+% subplot(5,4,7:8)
figure(1);
-subplot(5,4,7:8)
+subplot(3,4,5:6)%Revision 1
% subplot(2,4,7:8)
-for indexParam = 1:length(df2_global)
+snr_NAA_df2_global = zeros(length(df2_global),1);
+for indexParam = 1:length(df2_global)
for indexMetabolite = 1:numMetabolites
current_df2(indexMetabolite) = df2_global(indexParam);
end
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, water_ppm, scalingFactorRefPeak);
-
+
+ snr_NAA_df2_global(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFid))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, df2_global, titleStr, units);
end
% reset parameter
current_df2 = default_df2;
-titleStr = '$$\omega_{global}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-%% em - lorentzian components of T2
+%% em - lorentzian components of T2
paramKeyword = 'em';
+titleStr = '{\boldmath$\nu_{e}$}';
+units = '(Hz)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(df2_global), numMetabolites);
-figure(1);
-subplot(5,4,8+3:8+4)
+% figure(1);
+% subplot(5,4,8+3:8+4)
+figure(2);
+subplot(2,4,3:4)%Revision 1
% figure(2);
% subplot(3,4,3:4)
+snr_NAA_em = zeros(length(df2_global),1);
for indexParam = 1:length(df2_global)
for indexMetabolite = 1:numMetabolites
% generate a random number between [-stdLim; stdLim]
@@ -391,92 +424,136 @@ for indexParam = 1:length(df2_global)
case 'NAA(CH2)'
stdRand = NaaStdRand;
otherwise
- % nothing to do
+ % nothing to do
end
current_em_std_factor(indexMetabolite) = em_local_scale(indexParam) * stdRand;
end
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
-
+
+ snr_NAA_em(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFid))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, false, legendLabels, indexParam, em_local_scale, '$\rm T_{2,k} = T_{2,k} + X_{-1,1}^{k} \cdot $', '$\rm\cdot std(T_{2,k})$');
end
% reset parameter
current_em_std_factor = default_em_std_factor;
-titleStr = '{\boldmath$\nu_{e}$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
+legend(legendLabels, 'Interpreter', 'latex');
+plotArrowsPlotOffsets()
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-%% noise - add different noise levels
+%% noise - add different noise levels
paramKeyword = 'noise';
+titleStr = '{\boldmath$noise$}';
+units = '';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(noiseLevel), numMetabolites);
-figure(1);
-subplot(5,4,8+5:8+6)
+% figure(1);
+% subplot(5,4,8+5:8+6)
+figure(2);
+subplot(2,4,6:7)%Revision 1
% figure(2);
% subplot(3,4,5:6)
+snr_NAA_noise = zeros(length(noiseLevel),1);
for indexParam = 1:length(noiseLevel)
noiseVector = wgn(1,nTimepoints, noiseLevel(indexParam));
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
-
+
+ snr_NAA_noise(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, false, legendLabels, indexParam, snr_NAA_noise, '$\rm SNR_{NAA}$', units);
end
% reset parameter
noiseVector = wgn(1,nTimepoints, default_noiseLevel);
-titleStr = '{\boldmath$noise$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
+legend(legendLabels, 'Interpreter', 'latex');
+plotArrowsPlotOffsets()
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
-%% noise - add different noise levels
+%% noise - add different noise levels
paramKeyword = 'baseline';
+titleStr = '{\boldmath$baseline$}';
+units = '';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(noiseLevel), numMetabolites);
baselinesPath = [dataExportPathBase 'Baselines\'];
prefixBaseline = 'Baseline_';
baselines = {'1Flat_drop_H20'; '2Flat'; '3Flat'; '4Flat'; '5Wildish'; '6Wildish'; '7Very_wild';...
- '8Very_wild'; '9Wild'; '10Wild'; '11Wildish'; '12Wildish'; '13Normal'; '14Lipid_1.3'; '15Lipid_1.3_phased'; '16Zero'};
+ '8Very_wild'; '9Wild'; '10Wild'; '11Wildish'; '12Wildish'; '13Normal'; '14Lipid_1.3'; '15Lipid_1.3_phased'; '16Zero'};
+% figure(1);
+% subplot(5,4,8+10:8+11)
figure(1);
-subplot(5,4,8+10:8+11)
+subplot(3,4,11:12)
% figure(2);
% subplot(3,4,10:11)
+snr_NAA_baselines = zeros(length(baselines),1);
for indexParam = 1:length(baselines)
baselineFileName = strcat(prefixBaseline, baselines{indexParam});
[currentBaseline, ~, ~, ~, ~] = ImportLCModelBasis(baselinesPath, baselineFileName, plotBasis, '1H');
currentBaseline = currentBaseline(1:length(fids{1}));% will fail if the fids are longer than the baseline
- [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
+ [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak, currentBaseline);
-
+
+ snr_NAA_baselines(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
- hold on;
- plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
+
+ legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, baselines, titleStr, units);
end
-titleStr = '{\boldmath$baseline$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
end
+function legendLabels = plotSpectra(summedFidWithNoise, ppmVector, plotAll, legendLabels, indexParam, ...
+ paramValues, titleStr, units)
+
+if (plotAll == true)
+ %plot each fid
+ hold on;
+ plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
+else
+ if ~contains(titleStr,'=')
+ titleStr = [titleStr, ' = '];
+ end
+ % plot just the max, min and median value
+ offsetPlot = -1;
+ if (indexParam == 1)
+ offsetPlot = 0;
+ end
+ if (indexParam == round(length(paramValues)/2))
+ offsetPlot = 1;
+ end
+ if (indexParam == length(paramValues))
+ offsetPlot = 2;
+ end
+ if (offsetPlot ~= -1)
+ hold on;
+ plot(ppmVector,real(fftshift(fft(summedFidWithNoise)))-offsetPlot*2e3);
+ legendLabels{offsetPlot+1} = [titleStr, num2str(paramValues(indexParam)), ' ', units];
+ end
+end
+end
+
function plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase)
FontSize = 12;
LineWidth = 1.5;
@@ -498,3 +575,29 @@ if saveFigures
savefig(gcf, [fileNameFull, '.fig'],'compact');
end
end
+
+function plotArrowsPlotOffsets()
+
+arrowsColor = [0.4 0.4 0.4];
+offSetStep = 2*1e3;
+%arrows first - second
+yUpper = -0.2*offSetStep;
+yLower = -1.2*offSetStep;
+xa = [0.6 0.6]-0.05;
+doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
+%arrows second - third
+yUpper = -1.2*offSetStep;
+yLower = -2.2*offSetStep;
+xa = [0.6 0.6]-0.1;
+doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
+
+% text box with label "Plot offsets"
+xa = [0.6 0]+0.1;
+ya = [1e3 0];
+[xaf,yaf] = axescoord2figurecoord(xa,ya);
+str = {'Plot','offsets'};
+annotation('textbox',[xaf(1) yaf(1)-0.02 0.1 0.1],'String',str,'FitBoxToText','On', ...
+ 'EdgeColor','none', 'FontWeight','bold','Color',arrowsColor);
+
+set(gcf,'Color','w')
+end
diff --git a/createTestData/createSummedSpectrum.m b/createTestData/createSummedSpectrum.m
index 58bda11a3c7f8be6dbb15e9fb53c3a071840c137..a2be6d73fc47762389ebc1267050b85faed976dd 100644
--- a/createTestData/createSummedSpectrum.m
+++ b/createTestData/createSummedSpectrum.m
@@ -1,4 +1,4 @@
-function [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum( ...
+function [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum( ...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, shift_t2_mx, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyWord, indexParam, dwellTimeMs, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak, baseline)
@@ -13,7 +13,7 @@ currentConcentrationRm = size(1,length(metabolites));
for indexMetabolite = 1:length(metabolites)
current_T2 = T2s(indexMetabolite) + current_em_std_factor(indexMetabolite) * T2_std(indexMetabolite);%in ms
if strcmp(metabolites(indexMetabolite),'MMB')
- current_em(indexMetabolite) = 1;
+ current_em(indexMetabolite) = 0;
%the Gaussian for the acquired MMB should be smaller than the one for
%metabolites, since it was already measured with the inherent
%micro- and macro-susceptibility. It should however not be negative!
@@ -106,6 +106,9 @@ baselineScaled = baseline * maxValueSpectrum;
summedFidWithNoise = summedFid + noiseVector + baselineScaled;
% figure; plot(real(fftshift(fft(summedFid)))); hold on; plot(real(fftshift(fft(baseline*maxValueSpectrum))))
+bandwidth_Hz = 1 / dwellTimeMs * 1e3;
+stats = calculateSNR(scanFrequency, bandwidth_Hz, summedFidWithNoise);
+
if saveData
% export the spectra but also the used concentrations and all parameters
% export not truncated data
@@ -122,4 +125,89 @@ if saveData
save([dataExportPath, fileName, '.mat'], 'summedFidWithNoise', 'current_em', 'currentConcentrationRm', ...
'current_pc0', 'current_pc1', 'current_df2', 'current_gm', 'current_em_std_factor', 'conc_std_factor', 'baselineScaled');
end
+end
+
+function stats = calculateSNR(scanFrequency, bandwidth_Hz, fid)
+
+noiseRangePpm = [-4, -1];
+number_Pts = size(fid,2);
+
+bw = linspace(-bandwidth_Hz/2,bandwidth_Hz/2,number_Pts);
+ppmVector = bw/scanFrequency + 4.66;
+
+indexNoiseStart = find(ppmVector>noiseRangePpm(1),1,'first');
+indexNoiseEnd = find(ppmVector>noiseRangePpm(2),1,'first');
+noiseRange = [indexNoiseStart indexNoiseEnd];
+specReal = real(fftshift(fft(fid)));
+%calculate noise
+[noise_std] = noiseCalc(specReal, noiseRange);
+
+list = {[1.9 2.1] 'NAA'; [2.9 3.1] 'Cr';[2.25 2.5] 'Glu'; [3.1 3.3] 'tCho'};
+stats = cell(size(list,1)+1,1);
+stats{1,1} = 'metaboliteName';
+stats{1,2} = 'snr';
+stats{1,3} = 'linewidth (Hz)';
+stats{1,4} = 'amplitude';
+for i=1:size(list,1)
+ metRange = list{i,1};
+ metRange = convertPPMtoIndex(ppmVector, metRange);
+ metRange =[find(metRange==1,1,'first') find(metRange==1,1,'last')];
+
+ [snrV, amplitude] = snr(specReal,metRange,noise_std);
+ lineWidthValue = linewidth(specReal,metRange,bandwidth_Hz);
+ %area = sum(data(metRange(1):metRange(2)));
+ stats{i+1,1} = list{i,2};
+ stats{i+1,2} = snrV;
+ stats{i+1,3} = lineWidthValue;
+ stats{i+1,4} = amplitude;
+ %stats{i,4} = area;
+end
+end
+
+function [y] = convertPPMtoIndex(ppmVect, x)
+%CONVERTPPMTOINDEX converts the ppm pair x to indexes
+% Returns the vector of the size of the data points,
+% with the selected region set to 1
+ index1 = ppmVect >= x(1);
+ index2 = ppmVect <= x(2);
+ y = index1 & index2;
+end
+
+function [noise_std] = noiseCalc(fftsReal, noiseRange)
+% noise calculation
+noise = (fftsReal(noiseRange(1):noiseRange(2)))';
+fitobject = fit((0:length(noise)-1)',noise,'poly1'); %fit a poly2 to noise
+y = feval(fitobject,(0:length(noise)-1)'); % substract poly1 from noise
+noise = noise - y;% now the noise does not have linear component
+noise_std = std(noise);
+end
+
+function [snrValue, amplitude] = snr(specReal,peakRange,noise_std)
+amplitude = max(abs(specReal(peakRange(1):peakRange(2))));
+snrValue = amplitude/noise_std;
+%snrValue = 20*log(amplitude/noise); %in dB
+end
+
+function lineWidthValue = linewidth(specReal, metRange, bandwidth_Hz)
+peakRange = metRange;
+peak = specReal(peakRange(1):peakRange(2));
+
+interpFactor = 100; %
+indeces = [];
+i = 1;
+while (length(indeces)<2 && i <4)
+ peakTmp = interp(peak,interpFactor*i);
+ [~, maxPos] = max(abs(peakTmp));
+ normValues = peakTmp/(peakTmp(maxPos)/2);
+ indeces = find(normValues> 0.995 & normValues< 1.005);
+ i = i+1;
+end
+
+interpFactor = interpFactor*(i-1);
+freqResolution = bandwidth_Hz/(length(specReal)*interpFactor);
+try
+ lineWidthValue = (indeces(end) - indeces(1))*freqResolution;
+catch err
+ lineWidthValue = [];
+end
end
\ No newline at end of file
diff --git a/createTestData/doublePointedArrowInside.m b/createTestData/doublePointedArrowInside.m
new file mode 100644
index 0000000000000000000000000000000000000000..c28936fd065371c19e756230ff8e30bb0d6242a2
--- /dev/null
+++ b/createTestData/doublePointedArrowInside.m
@@ -0,0 +1,22 @@
+function doublePointedArrowInside(xa,yUpper,yLower, arrowsColor)
+ set(gcf,'Units','normalized');
+ yaDown = [yUpper+1e-8 yUpper];
+ yaLine = [yLower yUpper];
+ yaUp = [yLower-1e-8 yLower];
+ [xaf,yafLine] = axescoord2figurecoord(xa,yaLine);
+ line = annotation('line',xaf,yafLine);
+ [xaf,yafDown] = axescoord2figurecoord(xa,yaDown);
+ arDown = annotation('arrow',xaf,yafDown);
+ [xaf,yafUp] = axescoord2figurecoord(xa,yaUp);
+ arUp = annotation('arrow',xaf,yafUp);
+ %styling
+ line.Color = arrowsColor;
+ arDown.HeadStyle = 'vback2';
+ arDown.HeadLength = 5;
+ arDown.HeadWidth = 7;
+ arDown.Color = arrowsColor;
+ arUp.HeadStyle = 'vback2';
+ arUp.HeadLength = 5;
+ arUp.HeadWidth = 7;
+ arUp.Color = arrowsColor;
+end
\ No newline at end of file
diff --git a/createTestData/evaluateProFitFittingParameters.m b/createTestData/evaluateProFitFittingParameters.m
index df93d774d7475cf6d3e073d9bcece80449c8cb27..030808a7ffc2191591d2da4c73984e3601e1b136 100644
--- a/createTestData/evaluateProFitFittingParameters.m
+++ b/createTestData/evaluateProFitFittingParameters.m
@@ -4,8 +4,8 @@ folderName = 'final_simulations';
exportFolder = {'11 Simulations'};
pathNameExport = 'ProFit test data';
pathBaseExportFiles = pathToDataFolder(pathNameExport, exportFolder);
-% pathBaseExportFiles = pathBaseExportFiles(1:end-1);
-% pathBaseExportFiles = [pathBaseExportFiles ' - R3\'];
+pathBaseExportFiles = pathBaseExportFiles(1:end-1);
+pathBaseExportFiles = [pathBaseExportFiles ' - R3_new\'];% _only %+FID
dataExportPathBase = strcat(pathBaseExportFiles, exportFolder{1}, '\');
truncSuffix = '_truncOn';
folderFigures = [pathBaseExportFiles 'Output\' folderName, truncSuffix, '\'];
@@ -16,8 +16,8 @@ folders = {'df2', 'pc1', 'pc0', 'gauss'};
folders = {'df2_global', 'df2_local', 'pc1', 'pc0', 'gauss', 'em', 'noise', 'baseline'};
folders = {'pc0', 'pc1', 'gauss', 'df2_global', 'df2_local', ...
'em', 'noise', 'baseline'};
-xAxisLabels = {'$$\varphi_0$$', '$$\varphi_1$$', '$$\nu_g$$', '$$\omega_{global}$$', '{\boldmath$\omega_{local}$}', ...
- '{\boldmath$\nu_{e}$}', '{\boldmath$noise$}', '{\boldmath$baseline$}'};
+xAxisLabels = {'$$\rm\varphi_0$$', '$$\rm\varphi_1$$', '$$\rm\nu_g$$', '$$\rm\omega_{global}$$', '{\rm\boldmath$\omega_{local}$}', ...
+ '{\rm\boldmath$\nu_{e}$}', '{\rm\boldmath$noise$}', '{\rm\boldmath$baseline$}'};
legendLabels = {'\varphi_0 :\hspace{3em}', '\varphi_1 :\hspace{3.em}', '\nu_g :\hspace{3.em}', '\omega_{global} :\hspace{1.5em}', '${\boldmath$\omega_{local}$}$:\hspace{1.7em}', ...
'${\boldmath$\nu_{e}$}$:\hspace{3.5em}', '${\boldmath$noise$}$:\hspace{1.5em}', '${\boldmath$baseline$}$:'};
@@ -40,9 +40,15 @@ if useSubPlots
else
numParam = length(parameters);
end
+
+tableDeviations = cell(9,7);
+for iFolder = 1:length(folders)
+ tableDeviations{iFolder+1,1} = legendLabels{iFolder};
+end
for iParam = 1: numParam
legendLabels_ = {};
unitLabel = [titleLabels{iParam}, units{iParam}];
+ tableDeviations{1,iParam+1} = unitLabel;
for iFolder = 1:length(folders)
load([dataExportPathBase, folders{iFolder}, '\Diff_Summaries', truncSuffix, '_no_BL_df3.mat'], 'diff_conc_mat', 'diff_em_mat', 'diff_gm_mat', ...
'diff_pc0_mat', 'diff_df2_mat', 'diff_pc12_mat', 'diff_T2_mat');
@@ -71,19 +77,24 @@ for iParam = 1: numParam
else
scatter(zeros(1,size2)+iFolder, mean(diff_mat));
end
- title(['$$', titleLabels{iParam}, '$$'],'Interpreter','latex');
- legendLabels_{iFolder}= ['$$' legendLabels{iFolder}, num2str(diff_mean,'%.2f'), ' \pm ', num2str(diff_std,'%.2f'), '$$'];
- lgd = legend(legendLabels_,'Interpreter','latex','Location','northwest');
+ title(['$$\rm', titleLabels{iParam}, '$$'],'Interpreter','latex');
+ deviationsString = [ num2str(diff_mean,'%.2f'), ' \pm ', num2str(diff_std,'%.2f')];
+ tableDeviations{iFolder+1,iParam+1} = deviationsString;
+% legendLabels_{iFolder}= ['$$\rm' legendLabels{iFolder}, deviationsString, '$$'];
+% lgd = legend(legendLabels_,'Interpreter','latex','Location','northwest');
- title(lgd,{'$$Spectra\,\,with\,\,induced$$',['$$', titleOffset{iParam}, unitLabel,'$$']},'FontSize',8);
+% title(lgd,{'$$\rm Spectra\,\,with\,\,induced$$',['$$\rm', titleOffset{iParam}, unitLabel,'$$']},'FontSize',8);
end
xlim([0, length(folders)+1]);
xticks(1:length(folders));
xticklabels(xAxisLabels);
- xlabel('Simulations')
- set(gca,'YAxisLocation','right')
- ylabel(['$$', unitLabel, '$$'],'Interpreter','latex');
+% xlabel('Simulations')
+ xlabel('Spectra with induced variations in ')
+% set(gca,'YAxisLocation','right')
+ ylabel(['$$\rm', unitLabel, '$$'],'Interpreter','latex');
set(gca,'TickLabelInterpreter','latex')
set(gca,'XTickLabelRotation',60)
savefig([folderFigures, parameters{iParam}])
-end
\ No newline at end of file
+end
+
+writecell(tableDeviations,[folderFigures, 'Params.xlsx'])
\ No newline at end of file
diff --git a/createTestData/plotBaselines_and_fitted_spectra.m b/createTestData/plotBaselines_and_fitted_spectra.m
index 3d519c0652607eca535698b46254f40c1b20a026..100c6a07b480efbc714824961d9dca73feedd113 100644
--- a/createTestData/plotBaselines_and_fitted_spectra.m
+++ b/createTestData/plotBaselines_and_fitted_spectra.m
@@ -5,7 +5,7 @@ exportFolder = {'11 Simulations'};
pathNameExport = 'ProFit test data';
pathBaseExportFiles = pathToDataFolder(pathNameExport, exportFolder);
pathBaseExportFiles = pathBaseExportFiles(1:end-1);
-pathBaseExportFiles = [pathBaseExportFiles ' - R3\'];
+pathBaseExportFiles = [pathBaseExportFiles ' - R3_new\'];
dataExportPathBase = strcat(pathBaseExportFiles, exportFolder{1}, '\');
plotPath = [dataExportPathBase, 'Plots\BaselinePlots\'];
@@ -48,7 +48,8 @@ backgroundInlayColor = [0.8 0.8 0.8];
LineWidth = 1.5;
offsetBaseline = 7; % percent of max peak
-offsetResidual = 20; % percent of max peak
+offsetBaselineDiff = 18; % percent of max peak
+offsetResidual = 30; % percent of max peak
offsetPlots = 4; % percent of max peak
subplotGridX = 2;
subplotGridY = 2;
@@ -63,114 +64,125 @@ for indexParam = 1:length(baselines)
% figure;
figure('units','normalized','outerposition',[0.13,-0.014,0.85,0.99]); %Works on the Dell screen
end
- indexParamSorted = sortedFigure(indexParam);
- %% load data: Input data, LCModel fit, ProFit Fit
- % import the spectra but also the used concentrations and all parameters
- fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) '_truncOff'];
- load([dataExportPath, fileName, '.mat'], 'summedFidWithNoise', 'current_em', 'currentConcentrationRm', ...
- 'current_pc0', 'current_pc1', 'current_df2', 'current_gm', 'current_em_std_factor', 'conc_std_factor', 'baselineScaled');
- if plot_dkntmn
- fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) truncSuffix];
- [~, ~, phasedData_lcm_05, fitData_05, baselineData_lcm_05, residualData_05, ...
- ~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
- extractFits(dataExportPath_05,[fileName, '.coord']);
- [~, ~, phasedData_lcm_015, fitData_015, baselineData_lcm_015, residualData_015, ...
- ~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
- extractFits(dataExportPath_015,[fileName, '.coord']);
- end
- [nOfPoints, ppmVector_lcm, phasedData_lcm, fitData, baselineData_lcm, residualData, ...
- metaboliteNames, metaboliteData, tableConcentrations, ppmStart_lcm, ppmEnd_lcm, ...
- scanFrequency, frequencyShift, ppmGap, phase0_deg_lcm, phase1_deg_ppm_lcm] = ...
- extractFits(dataExportPath,[fileName, '.coord']);
-
- load(strcat(dataExportPath,fileName, '_profit.mat'),'fitresult', 'data');
-
- %% process input data
- nTimepoints = length(summedFidWithNoise);
- ppmVector = ppmVectorFunction(scanFrequency*1e6, bandwidth, nTimepoints, '1H');
+ indexParamSorted = sortedFigure(indexParam);
+ %% load data: Input data, LCModel fit, ProFit Fit
+ % import the spectra but also the used concentrations and all parameters
+ fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) '_truncOff'];
+ load([dataExportPath, fileName, '.mat'], 'summedFidWithNoise', 'current_em', 'currentConcentrationRm', ...
+ 'current_pc0', 'current_pc1', 'current_df2', 'current_gm', 'current_em_std_factor', 'conc_std_factor', 'baselineScaled');
+ if plot_dkntmn
+ fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) truncSuffix];
+ [~, ~, phasedData_lcm_05, fitData_05, baselineData_lcm_05, residualData_05, ...
+ ~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
+ extractFits(dataExportPath_05,[fileName, '.coord']);
+ [~, ~, phasedData_lcm_015, fitData_015, baselineData_lcm_015, residualData_015, ...
+ ~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
+ extractFits(dataExportPath_015,[fileName, '.coord']);
+ end
+ [nOfPoints, ppmVector_lcm, phasedData_lcm, fitData, baselineData_lcm, residualData, ...
+ metaboliteNames, metaboliteData, tableConcentrations, ppmStart_lcm, ppmEnd_lcm, ...
+ scanFrequency, frequencyShift, ppmGap, phase0_deg_lcm, phase1_deg_ppm_lcm] = ...
+ extractFits(dataExportPath,[fileName, '.coord']);
- spectrum = fftshift(fft(summedFidWithNoise));
- baselineSpectrum = fftshift(fft(baselineScaled));
-
- indexEnd = find(ppmVector > ppmEnd_lcm, 1);
- spectrumInputOffset = real(spectrum(indexEnd));
- baselineInputOffset = real(baselineSpectrum(indexEnd));
-
- %% process ProFit data
- currentFitresult = fitresult{1,iteration};
- ppm_profit = currentFitresult.common_data.data.fitroi.ppm2;
- %phasedData
- pData_profit = currentFitresult.processed_spec;
- %residual
- rData_profit = currentFitresult.residual_spec;
- %spline baseline
- if ~isempty(currentFitresult.splines_baseline)
- spBB = currentFitresult.splines_baseline.basis_matrix;
- concSpline = currentFitresult.splines_baseline.coeff;
- bData_profit = (spBB * concSpline(1:size(spBB,2)))';
- splines_lambda = currentFitresult.splines_lambda;
- splines_ed_vec = currentFitresult.splines_ed_vec;
- splines_optim_model_vec = currentFitresult.splines_optim_model_vec;
- splines_optim_idx = currentFitresult.splines_optim_idx;
- else
- bData_profit = zeros(size(currentFitresult.residual_spec));
- splines_lambda = 0;
- splines_ed_vec = [];
- splines_optim_model_vec = [];
- splines_optim_idx = 0;
- end
+ load(strcat(dataExportPath,fileName, '_profit.mat'),'fitresult', 'data');
+
+ %% process input data
+ nTimepoints = length(summedFidWithNoise);
+ ppmVector = ppmVectorFunction(scanFrequency*1e6, bandwidth, nTimepoints, '1H');
+
+ spectrum = fftshift(fft(summedFidWithNoise));
+ baselineSpectrum = fftshift(fft(baselineScaled));
+
+ indexEnd = find(ppmVector > ppmEnd_lcm, 1);
+ spectrumInputOffset = real(spectrum(indexEnd));
+ baselineInputOffset = real(baselineSpectrum(indexEnd));
+
+ %% process ProFit data
+ currentFitresult = fitresult{1,iteration};
+ ppm_profit = currentFitresult.common_data.data.fitroi.ppm2;
+ %phasedData
+ pData_profit = currentFitresult.processed_spec;
+ %residual
+ rData_profit = currentFitresult.residual_spec;
+ %spline baseline
+ if ~isempty(currentFitresult.splines_baseline)
+ spBB = currentFitresult.splines_baseline.basis_matrix;
+ concSpline = currentFitresult.splines_baseline.coeff;
+ bData_profit = (spBB * concSpline(1:size(spBB,2)))';
+ splines_lambda = currentFitresult.splines_lambda;
+ splines_ed_vec = currentFitresult.splines_ed_vec;
+ splines_optim_model_vec = currentFitresult.splines_optim_model_vec;
+ splines_optim_idx = currentFitresult.splines_optim_idx;
+ else
+ bData_profit = zeros(size(currentFitresult.residual_spec));
+ splines_lambda = 0;
+ splines_ed_vec = [];
+ splines_optim_model_vec = [];
+ splines_optim_idx = 0;
+ end
+
+ %% Scale for plotting
+ % calculate scaling. Maximum has to be from the same fitted range as in LCModel/ProFit
+ ppmMask = (ppmVector < ppmStart_lcm) & (ppmVector > ppmEnd_lcm);
+ maxSpectrum = max(real(spectrum(ppmMask)));
+ maxSpectrum_lcm = max(phasedData_lcm);
+ maxSpectrum_profit = max(real(pData_profit));
+ scalingFactor_lcm = maxSpectrum ./ maxSpectrum_lcm;
+ scalingFactor_profit = maxSpectrum ./ maxSpectrum_profit;
+
+ minBaseline = min(real(baselineSpectrum(ppmMask)));
+
+ %scale the offsets
+ offsetBaselineScaled = offsetBaseline * maxSpectrum / 100;
+ offsetBaselineDiffScaled = offsetBaselineDiff * maxSpectrum / 100 - minBaseline;
+ offsetResidualScaled = offsetResidual * maxSpectrum / 100 - minBaseline;
+ offsetPlotsScaled = offsetPlots * maxSpectrum / 100;
+
+ %% create the figure
+ indexSubplot = mod(indexParam-1, nGrid)+1;
+ hs = subplot(subplotGridX,subplotGridY,indexSubplot);
+ coordinates = hs.Position;
+
+ hold on; %if not alligned, mind that the ppmVector has to be defined with the water resonance at 4.7 ppm
+ plot(ppmVector,real(spectrum),'LineWidth', LineWidth, 'Color', colorSpectrum);
+ % plot(ppmVector_lcm,real(phasedData_lcm)*scalingFactor_lcm,'LineWidth', LineWidth, 'Color', colorLCModel);
+ % plot(ppm_profit,real(pData_profit)*scalingFactor_profit,'LineWidth', LineWidth, 'Color', colorProFit);
+
+ if plot_dkntmn
+ offsetBaselineScaled = offsetBaselineScaled + offsetPlotsScaled;
+ plot(ppmVector_lcm,real(baselineData_lcm_05)*scalingFactor_lcm-offsetBaselineScaled+2.4*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel_05);
+ plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+1.6*offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
+ plot(ppmVector_lcm,real(baselineData_lcm_015)*scalingFactor_lcm-offsetBaselineDiffScaled*2+0.8*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel);
+ else
+ plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
+ end
+ plot(ppmVector,real(baselineSpectrum)-offsetBaselineScaled,'-.','LineWidth', LineWidth, 'Color', colorOriginal);
+ plot(ppm_profit,real(bData_profit)*scalingFactor_profit-offsetBaselineScaled-offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorProFit);
+
+ fitBaselineLCM = fit(ppmVector_lcm, real(baselineData_lcm)*scalingFactor_lcm, 'linearinterp');
+ fitBaselineSpectrum = fit(ppmVector', real(baselineSpectrum)', 'linearinterp');
+ fitBaselineProFit = fit(ppm_profit', real(bData_profit)'*scalingFactor_profit, 'linearinterp');
+ plot(ppmVector_lcm, fitBaselineLCM(ppmVector_lcm)-fitBaselineSpectrum(ppmVector_lcm)-offsetBaselineDiffScaled,':','LineWidth', LineWidth, 'Color', colorLCModel)
+ plot(ppm_profit,fitBaselineProFit(ppm_profit)-fitBaselineSpectrum(ppm_profit)-offsetBaselineDiffScaled-offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorProFit);
+
+ plot(ppmVector_lcm,real(residualData)*scalingFactor_lcm-offsetResidualScaled,'LineWidth', LineWidth, 'Color', colorLCModel);
+ plot(ppm_profit,real(rData_profit)*scalingFactor_profit-offsetResidualScaled-offsetPlotsScaled,'LineWidth', LineWidth, 'Color', colorProFit);
+ set(gca,'xdir','reverse');
+ xlim([ppmEnd_lcm, ppmStart_lcm]);
+ ylim([-offsetResidualScaled-3*offsetPlotsScaled maxSpectrum+5*offsetPlotsScaled])
+ xlabel('\delta (ppm)');
+ ylabel('Signal (arb. u.)');
+ set(gca,'ytick',[]);
- %% Scale for plotting
- % calculate scaling. Maximum has to be from the same fitted range as in LCModel/ProFit
- ppmMask = (ppmVector < ppmStart_lcm) & (ppmVector > ppmEnd_lcm);
- maxSpectrum = max(real(spectrum(ppmMask)));
- maxSpectrum_lcm = max(phasedData_lcm);
- maxSpectrum_profit = max(real(pData_profit));
- scalingFactor_lcm = maxSpectrum ./ maxSpectrum_lcm;
- scalingFactor_profit = maxSpectrum ./ maxSpectrum_profit;
-
- minBaseline = min(real(baselineSpectrum(ppmMask)));
-
- %scale the offsets
- offsetBaselineScaled = offsetBaseline * maxSpectrum / 100;
- offsetResidualScaled = offsetResidual * maxSpectrum / 100 - minBaseline;
- offsetPlotsScaled = offsetPlots * maxSpectrum / 100;
-
- %% create the figure
- indexSubplot = mod(indexParam-1, nGrid)+1;
- hs = subplot(subplotGridX,subplotGridY,indexSubplot);
- coordinates = hs.Position;
-
- hold on; %if not alligned, mind that the ppmVector has to be defined with the water resonance at 4.7 ppm
- plot(ppmVector,real(spectrum),'LineWidth', LineWidth, 'Color', colorSpectrum);
-% plot(ppmVector_lcm,real(phasedData_lcm)*scalingFactor_lcm,'LineWidth', LineWidth, 'Color', colorLCModel);
-% plot(ppm_profit,real(pData_profit)*scalingFactor_profit,'LineWidth', LineWidth, 'Color', colorProFit);
-
- if plot_dkntmn
- offsetBaselineScaled = offsetBaselineScaled + offsetPlotsScaled;
- plot(ppmVector_lcm,real(baselineData_lcm_05)*scalingFactor_lcm-offsetBaselineScaled+2.4*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel_05);
- plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+1.6*offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
- plot(ppmVector_lcm,real(baselineData_lcm_015)*scalingFactor_lcm-offsetBaselineScaled+0.8*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel_015);
- else
- plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
- end
- plot(ppmVector,real(baselineSpectrum)-offsetBaselineScaled,'-.','LineWidth', LineWidth, 'Color', colorOriginal);
- plot(ppm_profit,real(bData_profit)*scalingFactor_profit-offsetBaselineScaled-offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorProFit);
-
- plot(ppmVector_lcm,real(residualData)*scalingFactor_lcm-offsetResidualScaled,'LineWidth', LineWidth, 'Color', colorLCModel);
- plot(ppm_profit,real(rData_profit)*scalingFactor_profit-offsetResidualScaled-offsetPlotsScaled,'LineWidth', LineWidth, 'Color', colorProFit);
- set(gca,'xdir','reverse');
- xlim([ppmEnd_lcm, ppmStart_lcm]);
- ylim([-offsetResidualScaled-3*offsetPlotsScaled maxSpectrum+5*offsetPlotsScaled])
- xlabel('\delta (ppm)');
- ylabel('Signal (arb. u.)');
- set(gca,'ytick',[]);
-
if plot_dkntmn
legend('Spectrum',sprintf('LCModel baseline\ndkntmn=0.5'),sprintf('LCModel baseline\ndkntmn=0.25'),sprintf('LCModel baseline\ndkntmn=0.15'),'Input baseline','ProFit baseline','LCModel residual','ProFit residual')
else
- legend('Spectrum','LCModel baseline','Input baseline','ProFit baseline','LCModel residual','ProFit residual')
+% legend('Spectrum','LCModel baseline','Input baseline','ProFit baseline','LCModel residual','ProFit residual')
+ if mod(indexParam,nGrid)==2
+ handle_legend = legend('Spectrum','LCModel baseline','Input baseline','ProFit-1D baseline',['LCModel ' char(8211) ' Input baseline'],['ProFit-1D ' char(8211) ' Input baseline'], 'LCModel residual','ProFit-1D residual','FontSize',12);
+ title(handle_legend,'Legend');
+ end
end
%arrows residual
@@ -188,7 +200,13 @@ for indexParam = 1:length(baselines)
%arrows Input baseline - ProFit
yUpper = baselineInputOffset-offsetBaselineScaled+offsetPlotsScaled;
yLower = baselineInputOffset-offsetBaselineScaled;
- xa = [ppmEnd_lcm ppmEnd_lcm]-0.1;
+ xa = [ppmEnd_lcm ppmEnd_lcm]-0.12;
+ doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
+
+ %arrows Baseline differences ProFit & LCModel
+ yUpper = baselineInputOffset-offsetBaselineDiffScaled;
+ yLower = baselineInputOffset-offsetBaselineDiffScaled-offsetPlotsScaled;
+ xa = [ppmEnd_lcm ppmEnd_lcm]-0.12;
doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
% text box with label "Plot offsets"
@@ -197,12 +215,13 @@ for indexParam = 1:length(baselines)
[xaf,yaf] = axescoord2figurecoord(xa,ya);
str = {'Plot','offsets'};
annotation('textbox',[xaf(1) yaf(1)-0.02 0.1 0.1],'String',str,'FitBoxToText','On', ...
- 'EdgeColor','none', 'FontWeight','bold','Color',arrowsColor);
+ 'EdgeColor','none', 'FontWeight','bold','Color',arrowsColor,'FontSize',13);
set(gcf,'Color','w')
+ set(gca,'FontSize',13)
- annotation('rectangle',[coordinates(1)+coordinates(3)*0.02 coordinates(2)+coordinates(4)*0.62 ...
- .35./subplotGridX .4./subplotGridY],'FaceColor', 'w', 'FaceAlpha',0, 'Color',[0 0.4 0.1])
+% annotation('rectangle',[coordinates(1)+coordinates(3)*0.02 coordinates(2)+coordinates(4)*0.62 ...
+% .35./subplotGridX .4./subplotGridY],'FaceColor', 'w', 'FaceAlpha',0, 'Color',[0 0.4 0.1])
axes('Position',[coordinates(1)+coordinates(3)*0.15 coordinates(2)+coordinates(4)*0.81 ...
.2./subplotGridX .2./subplotGridY])
@@ -210,20 +229,28 @@ for indexParam = 1:length(baselines)
plot(splines_ed_vec, splines_optim_model_vec,'Color',[0. 0.5 0.2]);
hold on
if ~isempty(splines_ed_vec(splines_optim_idx))
- xl = xline(splines_ed_vec(splines_optim_idx),'-.', {['ED: ' sprintf('%.f',splines_ed_vec(splines_optim_idx))],...
- ['mAIC: ', sprintf('%.1f',splines_optim_model_vec(splines_optim_idx))]}, ...
+% xl = xline(splines_ed_vec(splines_optim_idx),'-.', {['ED: ' sprintf('%.f',splines_ed_vec(splines_optim_idx))],...
+% ['mAIC: ', sprintf('%.1f',splines_optim_model_vec(splines_optim_idx))]}, ...
+% 'HandleVisibility','off', 'LineWidth', 2);
+% xl.LabelVerticalAlignment = 'top';
+% xl.LabelHorizontalAlignment = 'center';
+% xl.FontSize = 9;
+ xl = xline(splines_ed_vec(splines_optim_idx),'-.',...
'HandleVisibility','off', 'LineWidth', 2);
- xl.LabelVerticalAlignment = 'top';
- xl.LabelHorizontalAlignment = 'center';
- xl.FontSize = 9;
end
- set(gca,'FontSize',10)
+ set(gca,'FontSize',12)
set(gca,'Xscale','log');
% set(gca,'Color',backgroundInlayColor)
xlabel('Baseline ED / ppm');
ylabel('mAIC');
ylim([min(splines_optim_model_vec)*0.99 max(splines_optim_model_vec)*1.01])
- title('ProFit Baseline Smoothness');
+ % title('ProFit Baseline Smoothness');
+ if isempty(splines_ed_vec(splines_optim_idx))
+ title('ProFit Baseline Smoothness');
+ else
+ title(['ED: ' sprintf('%.f',splines_ed_vec(splines_optim_idx)),...
+ ' mAIC: ', sprintf('%.1f',splines_optim_model_vec(splines_optim_idx))]);
+ end
annotation('textbox',[coordinates(1)+coordinates(3)*0.02-0.03 coordinates(2)+coordinates(4)*0.62 ...
.35./subplotGridX .4./subplotGridY],'String',subplotNames{indexSubplot},'FontSize', 16, 'FontWeight', 'bold','LineStyle','none')
@@ -233,26 +260,3 @@ for indexParam = 1:length(baselines)
end
end
end
-
-function doublePointedArrowInside(xa,yUpper,yLower, arrowsColor)
- set(gcf,'Units','normalized');
- yaDown = [yUpper+1e-8 yUpper];
- yaLine = [yLower yUpper];
- yaUp = [yLower-1e-8 yLower];
- [xaf,yafLine] = axescoord2figurecoord(xa,yaLine);
- line = annotation('line',xaf,yafLine);
- [xaf,yafDown] = axescoord2figurecoord(xa,yaDown);
- arDown = annotation('arrow',xaf,yafDown);
- [xaf,yafUp] = axescoord2figurecoord(xa,yaUp);
- arUp = annotation('arrow',xaf,yafUp);
- %styling
- line.Color = arrowsColor;
- arDown.HeadStyle = 'vback2';
- arDown.HeadLength = 5;
- arDown.HeadWidth = 7;
- arDown.Color = arrowsColor;
- arUp.HeadStyle = 'vback2';
- arUp.HeadLength = 5;
- arUp.HeadWidth = 7;
- arUp.Color = arrowsColor;
-end
diff --git a/createTestData/plotSimulationResults.m b/createTestData/plotSimulationResults.m
index fde9ce3f5f8760969a910c63fd4c1f484c6c8de5..d461a5240715dab37363722a0af85c072e4f4d30 100644
--- a/createTestData/plotSimulationResults.m
+++ b/createTestData/plotSimulationResults.m
@@ -5,7 +5,7 @@ exportFolder = {'11 Simulations'};
pathNameExport = 'ProFit test data';
pathBaseExportFiles = pathToDataFolder(pathNameExport, exportFolder);
pathBaseExportFiles = pathBaseExportFiles(1:end-1);
-pathBaseExportFiles = [pathBaseExportFiles ' - R3\'];
+pathBaseExportFiles = [pathBaseExportFiles ' - R3_new\'];
dataExportPathBase = strcat(pathBaseExportFiles, exportFolder{1}, '\');
truncSuffix = '_truncOn';
scalingOn = true;
@@ -138,17 +138,19 @@ sumAllStdProFitMain = sumAllStdProFitMain + allStdProFitMain;
%% -----------------------------------------------------------------------------------------------------------
axes(figIds{1})
str = {'$$\varphi_0$$', '$$\varphi_1$$', '$$\omega_{global}$$', '{\boldmath$\omega_{local}$}'};
-xOffset = 0.65;
-yOffsetRect = 0.085; %0.075; -for two subfigures
+% xOffset = 0.65;
+xOffset = 0.38; %Revision 1
+yOffsetRect = 0.080; %0.085 - submission 1 %0.075; -for two subfigures
for indexSim = 1:4
- yOffset = 0.84;% 4 subfigs%0.80; - Two figs placed lower%0.86 -Two figs
+% yOffset = 0.84;% 4 subfigs%0.80; - Two figs placed lower%0.86 -Two figs
+ yOffset = 0.87;% Revision 1
annotation('textbox',[xOffset yOffset 0.1 0.1],'String',str{indexSim},'FitBoxToText','On', ...
'EdgeColor','none', 'FontSize', 12, 'Interpreter','latex','Color',plotColors{indexSim})
annotation('rectangle',[xOffset-0.012 yOffset+yOffsetRect 0.01 0.01],'EdgeColor','none', 'FaceColor',plotColors{indexSim})
- yOffset = 0.612;% 4 subfigs%0.37; - Two figs placed lower%0.37 -Two figs
- annotation('textbox',[xOffset yOffset 0.1 0.1],'String',str{indexSim},'FitBoxToText','On', ...
- 'EdgeColor','none', 'FontSize', 12, 'FontWeight','bold', 'Interpreter','latex','Color',plotColors{indexSim})
- annotation('rectangle',[xOffset-0.012 yOffset+yOffsetRect 0.01 0.01],'EdgeColor','none', 'FaceColor',plotColors{indexSim})
+% yOffset = 0.612;% 4 subfigs%0.37; - Two figs placed lower%0.37 -Two figs
+% annotation('textbox',[xOffset yOffset 0.1 0.1],'String',str{indexSim},'FitBoxToText','On', ...
+% 'EdgeColor','none', 'FontSize', 12, 'FontWeight','bold', 'Interpreter','latex','Color',plotColors{indexSim})
+% annotation('rectangle',[xOffset-0.012 yOffset+yOffsetRect 0.01 0.01],'EdgeColor','none', 'FaceColor',plotColors{indexSim})
xOffset = xOffset + 0.004*length(str{indexSim});
end
@@ -157,8 +159,9 @@ xlswrite([dataExportPathBase, 'SimulationsTable1.xlsx'],metaboliteTable);
numOfParameters2 = 5;
[metaboliteTable, indexTable] = setupTable(numberOfMet, numOfParameters2, metabolitesLabels);
-% figIds = [];
+figIds = [];
secondaryFigures = 2;
+secondaryFigures = 0;
%% -----------------------------------------------------------------------------------------------------------
%%
paramKeyword = 'gauss';
@@ -286,17 +289,19 @@ sumAllStdProFitMain = sumAllStdProFitMain / (numOfParameters + numOfParameters
axes(figIds{1})
str = {'$$\nu_g$$', '{\boldmath$\nu_{e}$}', '{\boldmath$noise$}', '{\boldmath$baseline$}'};
stringLenghts = {13, 13, 19, 27};
-xOffset = 0.63;
-yOffsetRect = 0.085; %0.075; -for two subfigures
+% xOffset = 0.63;
+xOffset = 0.38; %Revision 1
+yOffsetRect = 0.080; %0.085 - submission 1 %0.075; -for two subfigures
for indexSim = 1:4
- yOffset = 0.392;% 4 subfigs%0.80; - Two figs placed lower%0.86 -Two figs
+% yOffset = 0.392;% 4 subfigs%0.80; - Two figs placed lower%0.86 -Two figs
+ yOffset = 0.87;% Revision 1
annotation('textbox',[xOffset yOffset 0.1 0.1],'String',str{indexSim},'FitBoxToText','On', ...
'EdgeColor','none', 'FontSize', 12, 'Interpreter','latex','Color',plotColors{indexSim})
annotation('rectangle',[xOffset-0.012 yOffset+yOffsetRect 0.01 0.01],'EdgeColor','none', 'FaceColor',plotColors{indexSim})
- yOffset = 0.172;% 4 subfigs%0.37; - Two figs placed lower%0.37 -Two figs
- annotation('textbox',[xOffset yOffset 0.1 0.1],'String',str{indexSim},'FitBoxToText','On', ...
- 'EdgeColor','none', 'FontSize', 12, 'FontWeight','bold', 'Interpreter','latex','Color',plotColors{indexSim})
- annotation('rectangle',[xOffset-0.012 yOffset+yOffsetRect 0.01 0.01],'EdgeColor','none', 'FaceColor',plotColors{indexSim})
+% yOffset = 0.172;% 4 subfigs%0.37; - Two figs placed lower%0.37 -Two figs
+% annotation('textbox',[xOffset yOffset 0.1 0.1],'String',str{indexSim},'FitBoxToText','On', ...
+% 'EdgeColor','none', 'FontSize', 12, 'FontWeight','bold', 'Interpreter','latex','Color',plotColors{indexSim})
+% annotation('rectangle',[xOffset-0.012 yOffset+yOffsetRect 0.01 0.01],'EdgeColor','none', 'FaceColor',plotColors{indexSim})
xOffset = xOffset + 0.004*stringLenghts{indexSim};
end
@@ -341,6 +346,11 @@ allStdLCModelMain = mean(abs(stdConcDevLCModel(indexesMainMet)));
allMeanProFitMain = mean(abs(meanConcDevProFit(indexesMainMet)));
allStdProFitMain = mean(abs(stdConcDevProFit(indexesMainMet)));
+allMeanLCModelSmall = mean(abs(meanConcDevLCModel(~indexesMainMet)));
+allStdLCModelSmall = mean(abs(stdConcDevLCModel(~indexesMainMet)));
+allMeanProFitSmall = mean(abs(meanConcDevProFit(~indexesMainMet)));
+allStdProFitSmall = mean(abs(stdConcDevProFit(~indexesMainMet)));
+
plusMinusSign = char(177);
metaboliteTable{2,indexTable} = 'ProFit';
metaboliteTable{2,indexTable+1} = 'LCModel';
@@ -354,6 +364,9 @@ metaboliteTable{numberOfMet+3,indexTable+1} = [num2str(allMeanLCModel,'%.1f'),
metaboliteTable{numberOfMet+4,indexTable} = [num2str(allMeanProFitMain,'%.1f'), plusMinusSign, num2str(allStdProFitMain,'%.1f')];
metaboliteTable{numberOfMet+4,indexTable+1} = [num2str(allMeanLCModelMain,'%.1f'), plusMinusSign, num2str(allStdLCModelMain,'%.1f')];
+
+metaboliteTable{numberOfMet+5,indexTable} = [num2str(allMeanProFitSmall,'%.1f'), plusMinusSign, num2str(allStdProFitSmall,'%.1f')];
+metaboliteTable{numberOfMet+5,indexTable+1} = [num2str(allMeanLCModelSmall,'%.1f'), plusMinusSign, num2str(allStdLCModelSmall,'%.1f')];
end
%% reorganize the fitted concentrations to 2D tables
@@ -413,6 +426,12 @@ end
function figIdsCorr = plotCorrelationFigs(concentrationsRm, allLCModelMetConc, allProFitMetConc,...
metabolitesToDisplay, numberOfMet, dataExportPathBase, paramKeyword, truncSuffix, figIdsCorr)
+%TODO should I have this here? I don't know
+plotCorrelationsGlx(concentrationsRm, allLCModelMetConc, allProFitMetConc,...
+ metabolitesToDisplay,1:100);
+plotCorrelationsGlx(concentrationsRm, allLCModelMetConc, allProFitMetConc,...
+ metabolitesToDisplay,1:5:100);
+
if isempty(figIdsCorr)
emptyFigureIdsCorr = true;
figIdsCorr = cell(1,numberOfMet);
@@ -450,16 +469,23 @@ for indexMetabolite = 1:numberOfMet
figure(figID_Main);
mainFigureRunningIndex = mainFigureRunningIndex + 1;
subplot(figureMainRows, figureMainCols, mainFigureRunningIndex);
- if mod(mainFigureRunningIndex,figureMainCols)== 1
+% if mod(mainFigureRunningIndex,figureMainCols)== 1 %first submission. Label on every row
+ if mainFigureRunningIndex == 7 % Revision 1: Just one label
plotYLabel = true;
else
plotYLabel = false;
end
- if round((mainFigureRunningIndex+1)/figureMainCols)== figureMainRows
+% if round((mainFigureRunningIndex+1)/figureMainCols)== figureMainRows %first submission. Label on every column
+ if mainFigureRunningIndex == 11 %Revision 1: Just one label
plotXLabel = true;
else
plotXLabel = false;
end
+ if mainFigureRunningIndex == 3 %Revision 1: Just one legend
+ plotLegend = true;
+ else
+ plotLegend = false;
+ end
elseif sum(contains(metabolitesFigureSupporting, metabolitesToDisplay{indexMetabolite}))
figure(figID_Supporting);
supportingFigureRunningIndex = supportingFigureRunningIndex + 2;
@@ -477,6 +503,11 @@ for indexMetabolite = 1:numberOfMet
else
plotXLabel = false;
end
+ if supportingFigureRunningIndex == 4 %Revision 1: Just one legend
+ plotLegend = true;
+ else
+ plotLegend = false;
+ end
else
warning('Metabolite not associated with either main or supporting figure: ', metabolitesToDisplay{indexMetabolite});
figure();
@@ -498,13 +529,15 @@ for indexMetabolite = 1:numberOfMet
title(metabolitesToDisplay(indexMetabolite))
if plotXLabel
- xlabel('$c_{k}^{simulated} \,\, (mmol/kg)$', 'Interpreter', 'latex','FontSize', 10)
+ xlabel('$\bf c_{k}^{simulated} \,\, (mmol/kg)$', 'Interpreter', 'latex','FontSize', 11)
end
if plotYLabel
- ylabel('$c_k^{fitted} \,\, (mmol/kg)$', 'Interpreter', 'latex','FontSize', 10)
+ ylabel('$\bf c_k^{fitted} \,\, (mmol/kg)$', 'Interpreter', 'latex','FontSize', 11)
+ end
+ if plotLegend
+ handleLegend = legend('$LCModel \,\, c_k^{fitted}$', '$ProFit \,\, c_k^{fitted}$', '$Identity \, Line$', 'Location', 'best', 'Interpreter', 'latex', 'FontSize', 10);
+ title(handleLegend,'Legend', 'FontSize',10);
end
- legend('$LCModel \,\, c_k^{fitted}$', '$ProFit \,\, c_k^{fitted}$', '$Identity \, Line$', 'Location', 'best', 'Interpreter', 'latex', 'FontSize', 8)
-
xlim([min(identityLine),max(identityLine)])
ylim([min(identityLine),max(identityLine)])
@@ -525,47 +558,55 @@ function figIds = plotPrecisionFigs(figIds, numberOfMet, offsetPlot, secondaryFi
yLimValue = 100;
if isempty(figIds)
- % figure('Position',[612,509,1132,469],'OuterPosition',[604,501,1148,562], 'Units','pixels')
- % figure('Position',[100,275,890,680],'Units','pixels')
- % figIds{1} = subplot(2,1,2);
- % figIds{2} = subplot(2,1,1);%, 'Position', [0.58,0.16,0.35,0.76], 'Position', [0.13,0.16,0.35,0.76]
- %
- figure('Position',[-908,-589,890,1400],'Units','pixels')
- figIds{1} = subplot(4,1,2);
- figIds{2} = subplot(4,1,1);
- figIds{3} = subplot(4,1,4);
- figIds{4} = subplot(4,1,3);
+% figure('Position',[612,509,1132,469],'OuterPosition',[604,501,1148,562], 'Units','pixels')
+% figure('Position',[100,275,890,680],'Units','pixels')
+% figIds{1} = subplot(2,1,2);
+% figIds{2} = subplot(2,1,1);%, 'Position', [0.58,0.16,0.35,0.76], 'Position', [0.13,0.16,0.35,0.76]
+%%Revision 1
+% figure('Position',[612,509,1132,469],'OuterPosition',[604,501,1148,562], 'Units','pixels')
+ figure('Position',[9,10,1033,986],'Units','pixels')
+ figIds{1} = subplot(1,2,2);
+ figIds{2} = subplot(1,2,1);%, 'Position', [0.58,0.16,0.35,0.76], 'Position', [0.13,0.16,0.35,0.76]
+
+% figure('Position',[-908,-589,890,1400],'Units','pixels')
+% figIds{1} = subplot(4,1,2);
+% figIds{2} = subplot(4,1,1);
+% figIds{3} = subplot(4,1,4);
+% figIds{4} = subplot(4,1,3);
end
positionsTicks = [2.2:1.5:1.5*numberOfMet+1.5];
positions = positionsTicks + offsetPlot;
axes(figIds{1 + secondaryFigures})
hold on
- boxPlotLCModel = boxplot(allLCModelMetConcDev, metabolitesLabels,'colors',plotColors,'symbol','+', 'positions',positions,'width',0.18);
- xticks(positionsTicks)
- xticklabels(metabolitesLabels)
- xtickangle(45)
- xlim([1 27.5])
+ boxPlotLCModel = boxplot(allLCModelMetConcDev, metabolitesLabels,'colors',plotColors,'symbol','+', 'positions',positions,'width',0.18,'orientation', 'horizontal');
+ yticks(positionsTicks)
+ yticklabels(metabolitesLabels)
+ set(gca,'Ydir','reverse')
+% xtickangle(45)
title('LCModel')
- ylim([-yLimValue,yLimValue]);
- xlim([positionsTicks(1)-1, positionsTicks(end)+1])
- set(gca, 'FontSize', 11)
- ylabel('$Concentration\,\,change\,\,c^\%_{k}\,(\%)$', 'FontSize', 13, 'Interpreter', 'latex')
+ xlim([-yLimValue,yLimValue]);
+ ylim([positionsTicks(1)-1, positionsTicks(end)+1])
+ set(gca, 'FontSize', 12)
+ set(gca, 'YAxisLocation', 'right')
+ xlabel('$Concentration\,\,change\,\,c^\%_{k}\,(\%)$', 'FontSize', 13, 'Interpreter', 'latex')
set(boxPlotLCModel(:,:),'linewidth',1);
- yline(0,':', 'Color',[0.5 0.5 0.5]);
+ xline(0,':', 'Color',[0.5 0.5 0.5]);
%%
axes(figIds{2 + secondaryFigures})
hold on
- boxPlotProFit = boxplot(allProFitMetConcDev, metabolitesLabels,'colors',plotColors,'symbol','+', 'positions',positions,'width',0.18);
- xticks(positionsTicks)
- xticklabels(metabolitesLabels)
- xtickangle(45)
- xlim([1 27.5])
+ boxPlotProFit = boxplot(allProFitMetConcDev, metabolitesLabels,'colors',plotColors,'symbol','+', 'positions',positions,'width',0.18,'orientation', 'horizontal');
+ yticks(positionsTicks)
+ yticklabels(metabolitesLabels)
+ set(gca,'Ydir','reverse')
+% ytickangle(45)
+ ylim([1 27.5])
title('ProFit')
- ylim([-yLimValue,yLimValue]);
- set(gca, 'FontSize', 11)
- ylabel('$Concentration\,\,change\,\,c^\%_{k}\,(\%)$', 'FontSize', 13, 'Interpreter', 'latex')
+ xlim([-yLimValue,yLimValue]);
+ set(gca, 'FontSize', 12)
+ set(gca, 'YAxisLocation', 'left')
+ xlabel('$Concentration\,\,change\,\,c^\%_{k}\,(\%)$', 'FontSize', 13, 'Interpreter', 'latex')
set(boxPlotProFit(:,:),'linewidth',1);
- yline(0,':', 'Color',[0.5 0.5 0.5]);
+ xline(0,':', 'Color',[0.5 0.5 0.5]);
end
%% setup the table to show metabolite concentration changes
@@ -574,5 +615,62 @@ metaboliteTable = cell(numberOfMet+4,numOfParameters*2+1);
metaboliteTable(3:numberOfMet+2,1) = metabolitesLabels;
metaboliteTable{numberOfMet+3,1} = 'Mean';
metaboliteTable{numberOfMet+4,1} = 'Mean Main Metabolites';
+metaboliteTable{numberOfMet+5,1} = 'Mean Small Metabolites';
indexTable = 0;
+end
+
+function plotCorrelationsGlx(concentrationsRm, allLCModelMetConc, allProFitMetConc,...
+ metabolitesToDisplay, subsample)
+
+colorOriginal = [0 0 0];
+colorProFit = [0.49 0.18 0.56];
+colorLCModel = [0.85 0.33 0.1];
+
+indexGlu = find(strcmp(metabolitesToDisplay, 'Glu'));
+indexGln = find(strcmp(metabolitesToDisplay, 'Gln'));
+
+simulatedGluConc = concentrationsRm(:, indexGlu);
+simulatedGlnConc = concentrationsRm(:, indexGln);
+fittedGluConcLCModel = allLCModelMetConc(:, indexGlu);
+fittedGlnConcLCModel = allLCModelMetConc(:, indexGln);
+fittedGluConcProFit = allProFitMetConc(:, indexGlu);
+fittedGlnConcProFit = allProFitMetConc(:, indexGln);
+
+% subsample = 1:5:100;
+
+figure();
+plotYLabel = true;
+plotXLabel = true;
+plotLegend = true;
+
+scatterSize = 50;
+fontSize = 14;
+hold on
+scatter(fittedGluConcLCModel(subsample), fittedGlnConcLCModel(subsample), scatterSize, 'd', 'MarkerEdgeColor', colorLCModel);
+scatter(simulatedGluConc(subsample), simulatedGlnConc(subsample), scatterSize, 'o', 'MarkerEdgeColor', colorOriginal);
+scatter(fittedGluConcProFit(subsample), fittedGlnConcProFit(subsample), scatterSize, 'x', 'MarkerEdgeColor', colorProFit);
+for index = subsample
+ plot([fittedGluConcLCModel(index), simulatedGluConc(index), fittedGluConcProFit(index)],[fittedGlnConcLCModel(index), simulatedGlnConc(index), fittedGlnConcProFit(index)], 'color', colorOriginal);
+end
+
+offSetEnds = max(simulatedGluConc)*0.05;
+identityLineGlu = min(simulatedGluConc)-offSetEnds:0.01:max(simulatedGluConc)+offSetEnds;
+offSetEnds = max(simulatedGlnConc)*0.05;
+identityLineGln = min(simulatedGlnConc)-offSetEnds:0.01:max(simulatedGlnConc)+offSetEnds;
+
+set(gca, 'FontSize', fontSize)
+if plotXLabel
+ xlabel('$c_{Glu} \,\, (mmol/kg)$', 'Interpreter', 'latex','FontSize', fontSize+2)
+end
+if plotYLabel
+ ylabel('$c_{Gln} \,\, (mmol/kg)$', 'Interpreter', 'latex','FontSize', fontSize+2)
+end
+if plotLegend
+ handleLegend = legend('$LCModel$', '$Simulated$', '$ProFit$', 'Location', 'best', 'Interpreter', 'latex', 'FontSize', fontSize-2);
+ title(handleLegend,'Legend');
+end
+
+xlim([min(identityLineGlu),max(identityLineGlu)])
+ylim([min(identityLineGln),max(identityLineGln)])
+
end
\ No newline at end of file
diff --git a/createTestData/plot_all_test_data_profit.m b/createTestData/plot_all_test_data_profit.m
index bba246e127d4e0b6b907406286a063495378d13c..59431068a1a2343d3c913de4a1469dbeea444283 100644
--- a/createTestData/plot_all_test_data_profit.m
+++ b/createTestData/plot_all_test_data_profit.m
@@ -28,7 +28,7 @@ pathNameExport = 'ProFit test data';
pathBaseExportFiles = pathToDataFolder(pathNameExport, exportFolders);
pathBaseExportFiles = pathBaseExportFiles(1:end-1);
-pathBaseExportFiles = [pathBaseExportFiles ' - R3\'];
+pathBaseExportFiles = [pathBaseExportFiles ' - R3_new\'];
dataExportPathBaseRef = strcat(pathBaseExportFiles, exportFolders{2}, '\');
dataExportPathBase = strcat(pathBaseExportFiles, exportFolders{reconOption+1}, '\');
@@ -165,16 +165,17 @@ numOfSimulations = length(namingSuffixesAve)*length(namingSuffixesCoil);
% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'showFitCI',' on',...
% 'baStatsMode','Gaussian','forceZeroIntercept','on', 'diffValueMode', 'percent');
-rpc_table = cell(numberOfMet+5,5);
+rpc_table = cell(numberOfMet+6,5);
rpc_table{1,2} = softwareNames{1};
rpc_table{1,4} = softwareNames{2};
rpc_table{2,2} = subsets_refit{1};
rpc_table{2,3} = subsets_refit{2};
rpc_table{2,4} = subsets_refit{1};
rpc_table{2,5} = subsets_refit{2};
-rpc_table(3:end-3,1) = metabolitesToDisplay;
+rpc_table(3:end-4,1) = metabolitesToDisplay;
+rpc_table{end-2,1} = 'Mean';
rpc_table{end-1,1} = 'Mean Main';
-rpc_table{end,1} = 'Mean';
+rpc_table{end,1} = 'Mean Small';
%title('Reproducibility');
subplotsX = 4;
@@ -282,23 +283,39 @@ for indexMet = 1: numberOfMet
end
rpc_values = cell2mat(rpc_table(3:end-3,2:5));
+% extract summaries main met
rpc_meanProFit_32ave_main = mean(rpc_values(indecesMainMet,1));
rpc_meanProFit_64ave_main = mean(rpc_values(indecesMainMet,2));
rpc_meanLCModel_32ave_main = mean(rpc_values(indecesMainMet,3));
rpc_meanLCModel_64ave_main = mean(rpc_values(indecesMainMet,4));
-rpc_table(3:end-3,2:5) = sprintfc('%.0f',rpc_values);
-rpc_table{end-1,2} = num2str(rpc_meanProFit_32ave_main, '%.0f');
-rpc_table{end-1,3} = num2str(rpc_meanProFit_64ave_main, '%.0f');
-rpc_table{end-1,4} = num2str(rpc_meanLCModel_32ave_main, '%.0f');
-rpc_table{end-1,5} = num2str(rpc_meanLCModel_64ave_main, '%.0f');
+% extract summaries small met
+rpc_meanProFit_32ave_small = mean(rpc_values(~indecesMainMet,1));
+rpc_meanProFit_64ave_small = mean(rpc_values(~indecesMainMet,2));
+rpc_meanLCModel_32ave_small = mean(rpc_values(~indecesMainMet,3));
+rpc_meanLCModel_64ave_small = mean(rpc_values(~indecesMainMet,4));
+% extract summaries all met
rpc_meanProFit_32ave = mean(rpc_values(:,1));
rpc_meanProFit_64ave = mean(rpc_values(:,2));
rpc_meanLCModel_32ave = mean(rpc_values(:,3));
rpc_meanLCModel_64ave = mean(rpc_values(:,4));
-rpc_table{end,2} = num2str(rpc_meanProFit_32ave, '%.0f');
-rpc_table{end,3} = num2str(rpc_meanProFit_64ave, '%.0f');
-rpc_table{end,4} = num2str(rpc_meanLCModel_32ave, '%.0f');
-rpc_table{end,5} = num2str(rpc_meanLCModel_64ave, '%.0f');
+%% Add to RPC_table: all met
+rpc_table(3:end-4,2:5) = sprintfc('%.0f',rpc_values);
+% Add to RPC_table: mean values
+rpc_table{end-2,2} = num2str(rpc_meanProFit_32ave, '%.0f');
+rpc_table{end-2,3} = num2str(rpc_meanProFit_64ave, '%.0f');
+rpc_table{end-2,4} = num2str(rpc_meanLCModel_32ave, '%.0f');
+rpc_table{end-2,5} = num2str(rpc_meanLCModel_64ave, '%.0f');
+% Add to RPC_table: mean main met
+rpc_table{end-1,2} = num2str(rpc_meanProFit_32ave_main, '%.0f');
+rpc_table{end-1,3} = num2str(rpc_meanProFit_64ave_main, '%.0f');
+rpc_table{end-1,4} = num2str(rpc_meanLCModel_32ave_main, '%.0f');
+rpc_table{end-1,5} = num2str(rpc_meanLCModel_64ave_main, '%.0f');
+% Add to RPC_table: mean small met
+rpc_table{end,2} = num2str(rpc_meanProFit_32ave_small, '%.0f');
+rpc_table{end,3} = num2str(rpc_meanProFit_64ave_small, '%.0f');
+rpc_table{end,4} = num2str(rpc_meanLCModel_32ave_small, '%.0f');
+rpc_table{end,5} = num2str(rpc_meanLCModel_64ave_small, '%.0f');
+%write file
xlswrite([pathBaseExportFiles, 'InVivoResults_RPC.xlsx'], rpc_table)
%% code will break here. This used to be the old style plotting/results without Bland-Altman plots.
@@ -312,10 +329,15 @@ mainMetaboliteConcentrationsAllSubj = metaboliteConcentrationsAllSubjNorm(:, ind
meanMainMetDeviation = mean(abs(mainMetaboliteConcentrationsAllSubj(:)), 'omitnan');
stdMainMetDeviation = std(abs(mainMetaboliteConcentrationsAllSubj(:)), 'omitnan');
+smallMetaboliteConcentrationsAllSubj = metaboliteConcentrationsAllSubjNorm(:, ~indecesMainMet);
+meanSmallMetDeviation = mean(abs(smallMetaboliteConcentrationsAllSubj(:)), 'omitnan');
+stdSmallMetDeviation = std(abs(smallMetaboliteConcentrationsAllSubj(:)), 'omitnan');
+
metaboliteTable = cell(numberOfMet+3,2);
metaboliteTable(2:numberOfMet+1,1) = metabolitesLabels;
metaboliteTable{numberOfMet+2,1} = 'Mean';
metaboliteTable{numberOfMet+3,1} = 'Mean Main Metabolites';
+metaboliteTable{numberOfMet+4,1} = 'Mean Small Metabolites';
plusMinusSign = char(177);
metaboliteTable{1,2} = SoftwareName;
@@ -326,6 +348,7 @@ end
metaboliteTable{numberOfMet+2,2} = [num2str(meanAllMetDeviation,'%.1f'), plusMinusSign, num2str(stdAllMetDeviation,'%.1f')];
metaboliteTable{numberOfMet+3,2} = [num2str(meanMainMetDeviation,'%.1f'), plusMinusSign, num2str(stdMainMetDeviation,'%.1f')];
+metaboliteTable{numberOfMet+3,2} = [num2str(meanSmallMetDeviation,'%.1f'), plusMinusSign, num2str(stdSmallMetDeviation,'%.1f')];
xlswrite([pathBaseExportFiles, 'InVivoResults' SoftwareName, '.xlsx'], metaboliteTable)
%