Commit 77d0ccf6 authored by Antonio Politi's avatar Antonio Politi
Browse files

More error catching

parent c0dfc6fb
......@@ -34,3 +34,6 @@ Thumbs.db
matlabcode/FCS/example_data/optimisedFit/
focalVolume.txt
Focus.txt
matlabcode/FCS/example_data/Protein/
......@@ -287,7 +287,8 @@ classdef absFcsmodel < handle
tmp = C{1}(channelIdx(i));
[tokc] = regexp(tmp{:},'(\d+)', 'tokens');
%invert channel for compability with FA this assumes APD!!!
% invert channel for compability with FA this assumes APD!!!
% should check for ChS1 and ChS2
if str2num(char(tokc{:}))==2
ch = 1;
elseif str2num(char(tokc{:}))==1
......
General Correlation Corrections Correlation fitting
No File name Full path Format Clock Duration Channel mode Annotation Channel mode Selection 1st 2nd Frequency Eval frequency Mode Detrend frequency Offset Ch1 Amplitude Ch1 Rate Ch1 Offset Ch2 Amplitude Ch2 Rate Ch2 Xtalk Ch1 Xtalk Ch2 Labels Ch1 Labels Ch2 Interval1 1st Interval1 2nd Interval1 Ch1 Interval1 Ch2 Interval2 1st Interval2 2nd Interval2 Ch1 Interval2 Ch2 Bleaching Ch1 Bleaching Ch2 Bleaching Cross BG Ch1 BG Ch2 BG Cross Xtalk Ch1 Xtalk Ch2 Xtalk Cross Total Ch1 Total Ch2 Total Cross G0ave Ch1 G0std Ch1 G0ave Ch2 G0std Ch2 G0ave Cross G0std Cross CPMave Ch1 CPMstd Ch1 CPMave Ch2 CPMstd Ch2 CPMave Cross CPMstd Cross Model no Model name Mode 1st 2nd Chi sq Ch1 R sq adj Ch1 Chi sq Ch2 R sq adj Ch2 Chi sq Cross R sq adj Cross P01 name P01 value Ch1 P01 value Ch2 P01 value Cross P02 name P02 value Ch1 P02 value Ch2 P02 value Cross P03 name P03 value Ch1 P03 value Ch2 P03 value Cross P04 name P04 value Ch1 P04 value Ch2 P04 value Cross P05 name P05 value Ch1 P05 value Ch2 P05 value Cross P06 name P06 value Ch1 P06 value Ch2 P06 value Cross P07 name P07 value Ch1 P07 value Ch2 P07 value Cross P08 name P08 value Ch1 P08 value Ch2 P08 value Cross P09 name P09 value Ch1 P09 value Ch2 P09 value Cross P10 name P10 value Ch1 P10 value Ch2 P10 value Cross P11 name P11 value Ch1 P11 value Ch2 P11 value Cross P12 name P12 value Ch1 P12 value Ch2 P12 value Cross P13 name P13 value Ch1 P13 value Ch2 P13 value Cross P14 name P14 value Ch1 P14 value Ch2 P14 value Cross P15 name P15 value Ch1 P15 value Ch2 P15 value Cross P16 name P16 value Ch1 P16 value Ch2 P16 value Cross P17 name P17 value Ch1 P17 value Ch2 P17 value Cross P18 name P18 value Ch1 P18 value Ch2 P18 value Cross P19 name P19 value Ch1 P19 value Ch2 P19 value Cross P20 name P20 value Ch1 P20 value Ch2 P20 value Cross P21 name P21 value Ch1 P21 value Ch2 P21 value Cross P22 name P22 value Ch1 P22 value Ch2 P22 value Cross P23 name P23 value Ch1 P23 value Ch2 P23 value Cross P24 name P24 value Ch1 P24 value Ch2 P24 value Cross P25 name P25 value Ch1 P25 value Ch2 P25 value Cross
0 0001_R1_P1_K1_Ch2.zen C:\Users\toni\Dropbox\NPCMaturation\matlabcode\FCS\example_data\Protein\0001_R1_P1_K1_Ch2.zen FA/ISS 20000 30 1 NaN 1 0 0 29.94 100000 2000 Scale 0 1.72 0 0.01 0 0 0.01 0 0 1 1 0 1 488.16 0 0 1 488.16 0 1.0072 1 1 0.993 1 1 1 1 1 1.0001 1 1 0.0034939 0.00042527 0 0 0 0 1.705 0.20359 0 0 0 0 0 two-component anomalous diffusion with fluorescent protein-like blinking 1 0 0 0.0705007391736052 0.998487099311376 0 0 0 0 N 267.624847931722 10 10 thetaT 0.140406130654996 0.2 0.2 tauT 100 100 100 f1 1 1 1 tauD1 356.476754236024 500 500 alpha1 1.0030454268312 1 1 tauD2 5000 5000 5000 alpha2 1 1 1 kappa 3.5 5.5 5.5 offset 0 0 0 pValueLRT 0.271826311041273 0 0 pValueFtest 0.109209712107619 0 0 PAIC 1 0 0 BIC 29.5374984776539 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0
1 0001_R1_P2_K1_Ch2.zen C:\Users\toni\Dropbox\NPCMaturation\matlabcode\FCS\example_data\Protein\0001_R1_P2_K1_Ch2.zen FA/ISS 20000 30 1 NaN 1 0 0 29.94 100000 2000 Scale 0 1.72 0 0.01 0 0 0.01 0 0 1 1 0 1 495.126 0 0 1 495.126 0 0.9978 1 1 0.9931 1 1 1 1 1 0.9908 1 1 0.0038789 0.00076767 0 0 0 0 1.9172 0.35296 0 0 0 0 0 two-component anomalous diffusion with fluorescent protein-like blinking 1 0 0 0.0559507745430852 0.997093474324975 0 0 0 0 N 255.302757712356 10 10 thetaT 0.136358349521112 0.2 0.2 tauT 100 100 100 f1 1 1 1 tauD1 295.432609513387 500 500 alpha1 0.910428248032829 1 1 tauD2 5000 5000 5000 alpha2 1 1 1 kappa 3.5 5.5 5.5 offset 0 0 0 pValueLRT 0.729777263377255 0 0 pValueFtest 0.683955343795277 0 0 PAIC 1 0 0 BIC 26.6568995335047 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0
General Correlation Corrections Correlation fitting
No File name Full path Format Clock Duration Channel mode Annotation Channel mode Selection 1st 2nd Frequency Eval frequency Mode Detrend frequency Offset Ch1 Amplitude Ch1 Rate Ch1 Offset Ch2 Amplitude Ch2 Rate Ch2 Xtalk Ch1 Xtalk Ch2 Labels Ch1 Labels Ch2 Interval1 1st Interval1 2nd Interval1 Ch1 Interval1 Ch2 Interval2 1st Interval2 2nd Interval2 Ch1 Interval2 Ch2 Bleaching Ch1 Bleaching Ch2 Bleaching Cross BG Ch1 BG Ch2 BG Cross Xtalk Ch1 Xtalk Ch2 Xtalk Cross Total Ch1 Total Ch2 Total Cross G0ave Ch1 G0std Ch1 G0ave Ch2 G0std Ch2 G0ave Cross G0std Cross CPMave Ch1 CPMstd Ch1 CPMave Ch2 CPMstd Ch2 CPMave Cross CPMstd Cross Model no Model name Mode 1st 2nd Chi sq Ch1 R sq adj Ch1 Chi sq Ch2 R sq adj Ch2 Chi sq Cross R sq adj Cross P01 name P01 value Ch1 P01 value Ch2 P01 value Cross P02 name P02 value Ch1 P02 value Ch2 P02 value Cross P03 name P03 value Ch1 P03 value Ch2 P03 value Cross P04 name P04 value Ch1 P04 value Ch2 P04 value Cross P05 name P05 value Ch1 P05 value Ch2 P05 value Cross P06 name P06 value Ch1 P06 value Ch2 P06 value Cross P07 name P07 value Ch1 P07 value Ch2 P07 value Cross P08 name P08 value Ch1 P08 value Ch2 P08 value Cross P09 name P09 value Ch1 P09 value Ch2 P09 value Cross P10 name P10 value Ch1 P10 value Ch2 P10 value Cross P11 name P11 value Ch1 P11 value Ch2 P11 value Cross P12 name P12 value Ch1 P12 value Ch2 P12 value Cross P13 name P13 value Ch1 P13 value Ch2 P13 value Cross P14 name P14 value Ch1 P14 value Ch2 P14 value Cross P15 name P15 value Ch1 P15 value Ch2 P15 value Cross P16 name P16 value Ch1 P16 value Ch2 P16 value Cross P17 name P17 value Ch1 P17 value Ch2 P17 value Cross P18 name P18 value Ch1 P18 value Ch2 P18 value Cross P19 name P19 value Ch1 P19 value Ch2 P19 value Cross P20 name P20 value Ch1 P20 value Ch2 P20 value Cross P21 name P21 value Ch1 P21 value Ch2 P21 value Cross P22 name P22 value Ch1 P22 value Ch2 P22 value Cross P23 name P23 value Ch1 P23 value Ch2 P23 value Cross P24 name P24 value Ch1 P24 value Ch2 P24 value Cross P25 name P25 value Ch1 P25 value Ch2 P25 value Cross
0 0001_R1_P1_K1_Ch2.zen C:\Users\toni\Dropbox\NPCMaturation\matlabcode\FCS\example_data\Protein\0001_R1_P1_K1_Ch2.zen FA/ISS 20000 30 1 NaN 1 0 0 29.94 100000 2000 Scale 0 1.72 0 0.01 0 0 0.01 0 0 1 1 0 1 488.16 0 0 1 488.16 0 1.0072 1 1 0.993 1 1 1 1 1 1.0001 1 1 0.0034939 0.00042527 0 0 0 0 1.705 0.20359 0 0 0 0 0 two-component anomalous diffusion with fluorescent protein-like blinking 1 0 0 0.0424020359378806 0.999069241603751 0 0 0 0 N 273.770375377501 10 10 thetaT 0.187366006251353 0.2 0.2 tauT 100 100 100 f1 0.937210632693334 1 1 tauD1 363.159642143376 500 500 alpha1 1.19995735090635 1 1 tauD2 9659.22754143143 5000 5000 alpha2 1.15224595000963 1 1 kappa 3.5 5.5 5.5 offset 0 0 0 pValueLRT 0.271826311041273 0 0 pValueFtest 0.109209712107619 0 0 PAIC 0.26513042187774 0 0 BIC 40.4352005272802 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0
1 0001_R1_P2_K1_Ch2.zen C:\Users\toni\Dropbox\NPCMaturation\matlabcode\FCS\example_data\Protein\0001_R1_P2_K1_Ch2.zen FA/ISS 20000 30 1 NaN 1 0 0 29.94 100000 2000 Scale 0 1.72 0 0.01 0 0 0.01 0 0 1 1 0 1 495.126 0 0 1 495.126 0 0.9978 1 1 0.9931 1 1 1 1 1 0.9908 1 1 0.0038789 0.00076767 0 0 0 0 1.9172 0.35296 0 0 0 0 0 two-component anomalous diffusion with fluorescent protein-like blinking 1 0 0 0.0458942815275203 0.997556777884047 0 0 0 0 N 266.549942936583 10 10 thetaT 0.188998130240906 0.2 0.2 tauT 100 100 100 f1 0.944858003893859 1 1 tauD1 323.879316923508 500 500 alpha1 1.11984225696663 1 1 tauD2 10953.5688504368 5000 5000 alpha2 1.50531406954776 1 1 kappa 3.5 5.5 5.5 offset 0 0 0 pValueLRT 0.729777263377255 0 0 pValueFtest 0.683955343795277 0 0 PAIC 0.070269232138015 0 0 BIC 39.9390491475818 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0 NaN 0 0 0
......@@ -23,11 +23,15 @@ end
outputdir = 'optimisedFit';
% range to fit kappa
kappaVal = [2:0.1:7]';
kappaVal = [2:0.1:9]';
force = 0;
for idir = 1:length(indirCell)
indir = indirCell{idir};
display(['master_workflow_fit_dye processing directory ' num2str(idir) ' ' indir])
workflow_fit_dye(indir, kappaVal, outputdir, '1c');
try
workflow_fit_dye(indir, kappaVal, outputdir, '1c', force);
catch
display(['failed master_workflow_fit_dye processing directory ' num2str(idir) ' ' indir ' ']);
end
end
\ No newline at end of file
function master_workflow_fit_protein
function master_workflow_fit_protein(wd)
% MASTER_WORKFLOW_FIT_DYE call workflow_fit_dye for several directorie
% specified in the script
%
% Antonio Politi, EMBL, January 2017
wd = 'Z:\AntonioP_elltier1\CelllinesImaging\MitoSysPipelines';
force = 0;
%wd = 'Z:\AntonioP_elltier1\CelllinesImaging\MitoSysPipelines\160308_gfpNUP107z26z31\';
cd(wd);
if nargin <1
wd = uigetdir('.', 'Specify main directory where to search protein folders');
%wd = 'Z:\AntonioP_elltier1\CelllinesImaging\MitoSysPipelines';
cd(wd);
end
if ~exist(fullfile(wd, 'ProcessingHelpFiles'));
mkdir(fullfile(wd, 'ProcessingHelpFiles'));
end
......@@ -19,13 +20,21 @@ if ~exist(fullfile(wd, 'ProcessingHelpFiles', 'protein2cfiles.mat')) || force
for i = 1:length(indirCell)
resfiles = [resfiles; getAllFiles(indirCell{i}, 'res', '2c.res', 0)];
end
if isempty(resfiles)
resftable = readtable(fullfile(wd, 'ProcessingHelpFiles', 'protein2cfiles.txt'), 'ReadVariableName', 0);
for i=1:size(resftable,1)
resfiles = [resfiles;resftable(i,:).Var1];
end
resfiles = cell(resfiles);
end
save(fullfile(wd, 'ProcessingHelpFiles', 'protein2cfiles.mat'), 'resfiles');
else
load(fullfile(wd, 'ProcessingHelpFiles', 'protein2cfiles.mat'));
end
for ifile = 38:length(resfiles)
%%
%6 9 11
for ifile = 1:length(resfiles)
resfile = resfiles{ifile};
display(sprintf('master_workflow_fit_proteins processing res file %d/%d %s', ifile, length(resfiles), resfiles{ifile}))
fnames = workflow_fit_protein(resfile, '2c');
......@@ -93,7 +102,7 @@ function distributefiles(resfile, patterns)
tmpfile = fullfile(mdir{ip}, [ 'tmp.res']);
outfile = fullfile(mdir{ip}, [ basename '.res']);
writetable(res(idxs{ip},:), tmpfile , 'WriteVariableNames',false, 'Delimiter', '\t','FileType', 'text');
outfid= fopen(outfile, 'w');
outfid = fopen(outfile, 'w');
fprintf(outfid, header);
fclose(outfid);
appendFiles(tmpfile, outfile);
......
function [focVol] = workflow_fit_dye(indir, kappaVal, outputdir, FAsession)
function [focVol] = workflow_fit_dye(indir, kappaVal, outputdir, FAsession, force)
% COMPUTEDIRECTORY perform fitting of dye in a specific directory. program
% assumes that traces of one experiment are all in one directory!
% INPUT:
......@@ -31,7 +31,7 @@ if nargin < 1
end
if nargin < 2
kappaVal = [2:0.1:7]';
kappaVal = [2:0.1:9]';
end
if nargin < 3
......@@ -39,9 +39,12 @@ if nargin < 3
end
if nargin < 4
FAsession = '.1c.';
FAsession = '1c';
end
if nargin < 5
force = 0;
end
%% Initialization
% load model class
clear('MO');
......@@ -54,7 +57,7 @@ if ~exist(outdir)
end
% find files from ZEN and FA
zenfiles = getAllFiles(indir, 'fcs');
xmlFAfiles = getAllFiles(indir, 'xml', FAsession);
xmlFAfiles = getAllFiles(indir, 'xml',['\.' FAsession '\.']);
%corFiles = getAllFiles(indir,'cor', '1c'); this is obsolete as we do extract data from xml
%%
......@@ -62,7 +65,7 @@ if isempty(xmlFAfiles)
display('No FA xml found use .fcs file instead');
cond = 3;
else
cond = 1;
cond = 2;
end
switch cond
......@@ -97,13 +100,26 @@ end
norm = [];
%% Workflow
kappaValIn = kappaVal;
lk = length(kappaVal);
process = 1;
if exist(fullfile(outdir, [prefix '.mat'])) %load result file if it exists already
load(fullfile(outdir, [prefix '.mat']));
else
process = 0;
if length(Chfits(1).kappaVal) ~= length(kappaValIn)
display(sprintf('---------\n%s\n kappa boundaries have been changed. Data will be reprocessed.', indir));
process =1;
elseif (Chfits(1).kappaVal ~= kappaValIn)
display(sprintf('---------\n%s\n kappa boundaries have been changed. Data will be reprocessed.', indir));
process =1;
end
end
if (process || force)
% starting value of kappa
kappa = 5;
% number of repeats for each fit (starts from random values)
%MO.maxfittry = 1;
%maxRounds = 1;
MO.maxfittry = 5;
maxRounds = 3;
% struct to store data for each channel
......@@ -119,7 +135,7 @@ else
for irounds = 1:maxRounds % repeat workflow 3 times with updated optimal kappavalue at each round
%% find optimal tauT and thetaT
display(sprintf('Alexa dye fit optimal tauT and thetaT for Channel %d round %d/3', iC, irounds));
display(sprintf('Alexa dye fit optimal tauT and thetaT for Channel %d round %d/%d', iC, irounds, maxRounds));
%default starting values
MO.par.tauT = 7;
MO.par.thetaT = 0.2;
......@@ -129,7 +145,7 @@ else
MO.par.thetaT = thetaT;
%% fit N, tauD1, thetaT for constant tauT and different kappa values
display(sprintf('Alexa dye fit different kappa values for Channel %d round %d/3', iC, irounds));
display(sprintf('Alexa dye fit different kappa values for Channel %d round %d/%d', iC, irounds, maxRounds));
% [N_kappaVal*Nexp Npar] array of parameters (all) for each kappa
% value and experiment. Last value is norm
outPar = [];
......@@ -193,7 +209,8 @@ else
end
% avoid cases were boundary of kappaVal is touched
% avoid cases were boundary of kappaVal is touched and use a value of kappa
% = 5 instead
if all([Chfits.posBest] == [1 1]) || all([Chfits.posBest] == [lk lk])
posB = find(kappaVal==5);
......@@ -225,11 +242,12 @@ for iC = 1:2
fprintf(fid,'N\tthetaT\ttauT\tf1\ttauD1\talpha1\ttauD2\talpha2\tkappa\tnorm\tR2\tNpt\tid_exp\r\n');
outPar = Chfits(iC).outPar;
for i = 1:size(outPar,1)
for j = 1:size(outPar,2)
fprintf(fid,'%.4f\t', outPar(i,j));
for j = 1:size(outPar,2)-1
fprintf(fid,'%.4e\t', outPar(i,j));
end
fprintf(fid,'\r\n');
fprintf(fid,'%.4e\r\n', outPar(i,end));
end
fclose(fid);
end
%% update xml files and write a res file
......@@ -243,7 +261,7 @@ if ~ZEN
fit = zeros(Nfits, 3);
par = repmat(struct2array(MO.par),3,1);
nrPar = length(par);
for iC = 1:3
for iC = 1:2
if (dataC{iD}(1,iC*2) == 0) % do not process if channel has not been acquired
continue
end
......
......@@ -10,21 +10,21 @@ function outfiles = workflow_fit_protein(resfile, FAsession)
%% default values
if nargin < 1
% prompt for a directory
[resfile, pathname] = uigetfile({'*.res'; '*.*'},'Select res file', 'C:\Users\toni\Dropbox\NPCMaturation\matlabcode\FCS\example_data\Protein');
[resfile, pathname] = uigetfile({'*.res'; '*.*'},'Select res file', '');
%[resfile, pathname] = uigetfile({'*.res'; '*.*'},'Select res file', 'Z:\AntonioP_elltier1\CelllinesImaging\MitoSysPipelines\160308_gfpNUP107z26z31\MitoSys2');
if ~resfile
return
end
resfile = fullfile(pathname, resfile);
end
[pathname, resbase, ext] = fileparts(resfile);
flintfile = fullfile(pathname, [resbase '.fint']);
if nargin < 2
FAsession = '2c';
end
FAsession = resbase;
%FAsession = '2c'
end
if ~exist(resfile)
display(['Could not find resfile' resfile]);
......@@ -44,16 +44,21 @@ end
restable2c = readtable(resfile, 'Delimiter', '\t', 'Header', 1, 'FileType', 'text', 'TreatAsEmpty',{'NA'});
restable1c = restable2c;
% stat entries
restable1c.P11Name = repmat({'pValueLRT'}, length(restable1c.P11Name),1);
restable1c.P12Name = repmat({'pValueFtest'}, length(restable1c.P11Name),1);
restable1c.P13Name = repmat({'PAIC'}, length(restable1c.P11Name),1);
restable1c.P14Name = repmat({'BIC'}, length(restable1c.P11Name),1);
% focal volume entries
restable1c.P11Name = repmat({'w0'}, length(restable1c.P11Name),1);
restable1c.P12Name = repmat({'Veff'}, length(restable1c.P11Name),1);
restable2c.P11Name = repmat({'w0'}, length(restable1c.P11Name),1);
restable2c.P12Name = repmat({'Veff'}, length(restable1c.P11Name),1);
restable2c.P11Name = repmat({'pValueLRT'}, length(restable1c.P11Name),1);
restable2c.P12Name = repmat({'pValueFtest'}, length(restable1c.P11Name),1);
restable2c.P13Name = repmat({'PAIC'}, length(restable1c.P11Name),1);
restable2c.P14Name = repmat({'BIC'}, length(restable1c.P11Name),1);
% stat entries to compare to 2c to 1c
restable1c.P13Name = repmat({'pValueLRT'}, length(restable1c.P11Name),1);
restable1c.P14Name = repmat({'pValueFtest'}, length(restable1c.P11Name),1);
restable1c.P15Name = repmat({'PAIC'}, length(restable1c.P11Name),1);
restable1c.P16Name = repmat({'BIC'}, length(restable1c.P11Name),1);
restable2c.P13Name = repmat({'pValueLRT'}, length(restable1c.P11Name),1);
restable2c.P14Name = repmat({'pValueFtest'}, length(restable1c.P11Name),1);
restable2c.P15Name = repmat({'PAIC'}, length(restable1c.P11Name),1);
restable2c.P16Name = repmat({'BIC'}, length(restable1c.P11Name),1);
%% write fluorescence intensity to table
if ~isempty(flintT)
restable1c.class = flintT.class;
......@@ -71,20 +76,27 @@ end
%%
% get names of xml files
xmlFAfiles = strrep(restable2c.FullPath, ['.zen'], ['.zen.' FAsession '.xml']);
%xmlFAfiles = strrep(restable2c.FullPath, ['.zen_mod.fcs'], ['.zen_mod.fcs.' FAsession '.xml']);
% name should be with respect to subfolders! TODO
%% read focal volume data
focVol = getAllFiles(fileparts(resfile), 'txt', 'focalVolume');
if length(focVol) == 0
display(['No focalVolume file associated with ' resfile]);
return;
[focVol, pathname] = uigetfile({'*.txt'; '*.*'},'Select focal volume file', fileparts(resfile));
if ~focVol
return
end
focVol = {fullfile(pathname,focVol)};
end
if length(focVol) > 1
display(['several focalVolume files associated with ' resfile]);
return;
end
focVol = readtable(focVol{1}, 'Delimiter', '\t');
%%
if isempty(xmlFAfiles)
display('No FA xml found use .fcs file instead');
......@@ -118,13 +130,14 @@ end
MO.maxfittry = 3; %3 repetitions for each fit
MO.weightChi2 = 0;
weightStr = {'','_w'};
% in FA seems more that unweighted fit is used. The residuals look most
% in FA seems that unweighted fit is used. The residuals between matlab fit and FA fit look
% close to each other when unweighted is used. Also the first point is
% often out of fit. Unfortunately the exact measure on how the fit is
% performed is not available. If we take best fit here and plugin into FA
% the Chi2 slightly increases
for iD = 1:length(dataC)
display(sprintf('Fitting of data set %d/%d', iD, length(dataC)));
% Fit for all channels (if CC or XC is found)
......@@ -137,9 +150,13 @@ for iD = 1:length(dataC)
else
Chstr = ['Ch' num2str(iCh)];
end
if isempty(dataC{iD})
continue
end
if dataC{iD}(1,iCh*2) == 0 % exclude if no data available
continue
end
MO.par.kappa = focVol.kappa(iCh);
% norm1c is sum of squared residuals used in the fit. norms1c = [unweighted, weighted norm]
% R2_1c is from norm used for the fit. R2s_1c:[R2_unweighted,
......@@ -166,20 +183,28 @@ for iD = 1:length(dataC)
%report unweighted R2 as the weight are largely overestimated. This
restable1c(iD, ['RSqAdj' Chstr]) = table(R2adj_1c(1));
restable2c(iD, ['RSqAdj' Chstr]) = table(R2adj_2c(1));
restable1c(iD, ['P11Value' Chstr]) = table(focVol.w0_um(iCh));
restable1c(iD, ['P12Value' Chstr]) = table(focVol.Vol_fl(iCh));
restable2c(iD, ['P11Value' Chstr]) = table(focVol.w0_um(iCh));
restable2c(iD, ['P12Value' Chstr]) = table(focVol.Vol_fl(iCh));
restable1c(iD, ['P11Value' Chstr]) = table(pvalueChi);
restable1c(iD, ['P12Value' Chstr]) = table(pvalueFtest);
restable1c(iD, ['P13Value' Chstr]) = table(PAIC(1));
restable1c(iD, ['P14Value' Chstr]) = table(BIC(1));
restable2c(iD, ['P11Value' Chstr]) = table(pvalueChi);
restable2c(iD, ['P12Value' Chstr]) = table(pvalueFtest);
restable2c(iD, ['P13Value' Chstr]) = table(PAIC(2));
restable2c(iD, ['P14Value' Chstr]) = table(BIC(2));
restable1c(iD, ['P13Value' Chstr]) = table(pvalueChi);
restable1c(iD, ['P14Value' Chstr]) = table(pvalueFtest);
restable1c(iD, ['P15Value' Chstr]) = table(PAIC(1));
restable1c(iD, ['P16Value' Chstr]) = table(BIC(1));
restable2c(iD, ['P13Value' Chstr]) = table(pvalueChi);
restable2c(iD, ['P14Value' Chstr]) = table(pvalueFtest);
restable2c(iD, ['P15Value' Chstr]) = table(PAIC(2));
restable2c(iD, ['P16Value' Chstr]) = table(BIC(2));
end
% write updated xml file
for ifit = 1:2
if isempty(dataC{iD})
continue
end
node = xmlFA.readnode(xmlFAfiles{iD});
[path1, fout, ext] = fileparts(xmlFAfiles{iD});
[path, fout, ext] = fileparts(fout);
......@@ -209,10 +234,10 @@ end
%%
tmpfiles = {fullfile(fileparts(resfile), ['1c_tmp.res']), ...
fullfile(fileparts(resfile), ['2c_tmp.res'])};
fullfile(fileparts(resfile), [FAsession '_tmp.res'])};
outfiles = {fullfile(fileparts(resfile), ['1c_opt' weightStr{MO.weightChi2+1} '.res']), ...
fullfile(fileparts(resfile), ['2c_opt' weightStr{MO.weightChi2+1} '.res'])};
fullfile(fileparts(resfile), [FAsession '_opt' weightStr{MO.weightChi2+1} '.res'])};
tables = {restable1c, restable2c};
for ifile = 1:2
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment