Skip to content

GitLab

Commit 9ac08059 authored by Antonio Politi's avatar Antonio Politi
Browse files

modify indexing of boundary

parent 5ec061d4
Hide whitespace changes
Inline Side-by-side
matlabcode/FCS/absFcsmodel.m
View file @ 9ac08059
......@@ -16,6 +16,8 @@ classdef absFcsmodel < handle
par; % parameter values a struct with name of parameter and value. This needs to be defined in the derived class
lb; % low boundary of parameter values
hb; % high boundary of parameter values
tauBoundary;
model;
end
......@@ -62,7 +64,30 @@ classdef absFcsmodel < handle
end
end
function [N, tauD1, tauD2] = initialGuess(MO, data, tauStart, tauEnd)
function [idxs] = getDataIndexes(MO, data, tauBoundary)
% GETDATAINDEXES get indexes within tauBoundary
% data - is [tau, corr] array
% tauBoundary - is [], [tauStart], [tauStart, tauEnd]
% idxs - are indexes of where to start (base 1) and
% where to end. idxs = [1, 0] use whole data set. Compatible
% with FA notation
tauStart = data(1,1);
tauEnd = data(end,1);
if nargin == 3
if ~isempty(tauBoundary)
tauStart = tauBoundary(1);
if length(tauBoundary) == 2
tauEnd = tauBoundary(2);
end
end
end
ist= find(data(:,1) >= tauStart); % index of start time
ien = find(data(:,1) >= tauEnd); % index of end time
idxs(1) = ist(1);
idxs(2) = size(data,1) - ien(1);
end
function [N, tauD1, tauD2] = initialGuess(MO, data, tauBoundary)
% INITIALGUESS given the AC data returns an initial guess forN, tauD1 of first component, tauD2 of
% second component
% [N, tauD1, tauD2] = MO.initialGuess(data, tauStart, tauEnd)
......@@ -72,14 +97,21 @@ classdef absFcsmodel < handle
% Method from Wachsmuth, M., et al. (2015). Nature Biotechnology, 33(4), 384?389.
% https://doi.org/10.1038/nbt.3146
ist= find(data(:,1) >= tauStart); % index of start time
ien = find(data(:,1) >= tauEnd); % index of end time
ist = ist(1);
ien = ien(1);
N = 1/mean(data(ist:ist+5,2)); % get initial guess from average of AC first 5 time pts
N = N + (1-2*rand)*N/8; % randomize initial guess
if nargin < 3
tauBoundary = MO.tauBoundary;
end
idxboundary = MO.getDataIndexes(data, tauBoundary);
data = data(idxboundary(1):end-idxboundary(2),:);
% only positive values for N and tauDm
idx = find(data(:,2)>0);
data = data(idx, :);
N = 1/mean(data(1:5,2));
% compute integral of equation to obtain a gues of average time
% only positive values
tauDm = N*trapz(data(ist:ien,1),data(ist:ien,2))/4.5; % compute integral of equation to obtain a gues of average time
tauDm = N*trapz(data(:,1), data(:,2))/4.5;
%randomize
N = N + (1-2*rand)*N/8;
% another way of computing tauDm would be by using a linear approximation of AC
tauDm = tauDm + (1-2*rand)*tauDm/8; % randomize tauDm
% assign value of tauD1 and tauD2 so that tauD1 ~ 1/10 of tauD2. 1/tauDm = f1/tauD1 + (1-f1)/tauD2
......@@ -87,11 +119,13 @@ classdef absFcsmodel < handle
tauD2 = (9*MO.par.f1+1)*tauDm;
end
function Gdiff = Gdiffcomp(MO, par, data)
function Gdiff = Gdiffcomp(MO, par, data, weight)
% GDIFFCOMP computes weighted difference of ACF data and model
% par - full parameter struct of model data - Ntx3 matrix of data [time, ACF, ACF_std]
if MO.weightChi2
if nargin < 4
weight = MO.weightChi2;
end
if weight
Gdiff = (MO.Gcomp(par, data(:,1)) - data(:,2))./data(:,3);
else
Gdiff = (MO.Gcomp(par, data(:,1)) - data(:,2));
......@@ -110,18 +144,67 @@ classdef absFcsmodel < handle
G = MO.Gcomp(par, tau);
end
function Gdiff = Gdiffcompfit(MO, parf, names, par, data)
function Gdiff = Gdiffcompfit(MO, parf, names, par, data, weight)
% GDIFFCOMPFIT computes difference model to ACF using parameter specified for fit
% parf - parameter vector names - name of parameter to be changed par - default parameter
% struct tau - vector of ACF times
for i = 1:length(names)
par = setfield(par, names{i}, parf(i));
end
Gdiff = MO.Gdiffcomp(par, data);
if nargin < 6
weight = MO.weightChi2;
end
Gdiff = MO.Gdiffcomp(par, data, weight);
end
function [par, norm] = getbestfit(MO, parCA,normV)
% GETBESTFIT returns parameter and norm of best fit from
% several fits
% parCA - a cell array
% normV - a vector containing the norm
[val, pos] = min(normV);
norm = normV(pos);
par = parCA{pos};
end
function [R2, SSres] = compR2(MO, data, par, tauBoundary, weight)
%%
% COMPR2 computes R2 given the parameters
% Note that we are fitting a non-linear model so this has not
% much meaning for comparing models! It is just used to asses
% curves that are kind of flat.
% R2adj is used in FA
if nargin < 4
tauBoundary = MO.tauBoundary;
end
if nargin < 5
weight = MO.weightChi2;
end
idxboundary = MO.getDataIndexes(data, tauBoundary);
data = data(idxboundary(1):end-idxboundary(2),:);
if isstruct(par)
SSres = sum(MO.Gdiffcomp(par, data, weight).^2);
else
for i=1:size(par,1)
parS = MO.parvecToparstruct(par(i,:));
SSres(i,1) = sum(MO.Gdiffcomp(parS, data, weight).^2);
end
end
if weight
SStot = sum((data(:,2)./data(:,3) - 1./data(:,3)*mean(data(:,2))).^2);
else
SStot = sum((data(:,2) - mean(data(:,2))).^2);
end
R2 = 1- SSres./SStot;
% R2adj is corrected for the number of parameters that have
% been fitted
% R2adj_2c = 1 - (1-R2_2c)*(Nrpts2c-1)/(Nrpts2c - nrpara2c - 1) (https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2)
% We don't pass the Npar fitted to the function here
end
function parS = parvecToparstruct(MO, parV)
% PARVECTOPARSTRUCT return a struct of type MO.par with values given by parV
......@@ -140,13 +223,13 @@ classdef absFcsmodel < handle
semilogx(data(:,1), G,'-', data(:,1), data(:,2),'o');
end
function [par, norm] = fitmodel(MO, parf, names, par, data, tauStart, tauEnd)
% FITMODELfit model to data depending weightChi2 this use weighted difference or not
function [par, norm] = fitmodel(MO, parf, names, par, data, tauBoundary)
% FITMODEL fit model to data depending weightChi2 this use weighted difference or not
% [par, norm] = MO.fitmodel(parf, names, par, data) - fit in default time interval (set by the
% data) [
% [par, norm] = MO.fitmodel(parf, names, par, data, tauStart) - fit in time intervael
% [tauStart, tauMax] [par, norm] = MO.fitmodel(parf, names, par, data, tauStart, tauEnd) - fit in
% time intervael [tauStart, tauEnd]
% [par, norm] = MO.fitmodel(parf, names, par, data, tauBoundary) - fit in time intervael
% [tauBoundary(1), tauMax] or [tauBoundary(1) tauBoundary(2)]
% if tauBoundary has length 2
% INPUT:
% parf - vector of parameter values for the fit names - cell array with names of parameters.
% par - struct of parameters data - Ntx3 matrix of data [time, ACF, ACF_std] tauStart -
......@@ -154,16 +237,13 @@ classdef absFcsmodel < handle
% OUTPUT:
% par - struct containing fitted parameters norm - sum of squared difference at fitting
% minima
if nargin <6
tauStart = data(1,1);
end
if nargin < 7
tauEnd = data(end,1);
% norm - sum of squared residuals at minima
if nargin < 6
tauBoundary = MO.tauBoundary;
end
ist= find(data(:,1) >= tauStart);
ien = find(data(:,1) >= tauEnd);
ist = ist(1);
ien = ien(1);
idxboundary = MO.getDataIndexes(data, tauBoundary);
data = data(idxboundary(1):end-idxboundary(2),:);
lb = [];
hb = [];
for i = 1:length(names)
......@@ -171,7 +251,7 @@ classdef absFcsmodel < handle
hb = [hb;getfield(MO.hb,names{i})];
end
options = optimset('display','off');
[parf, norm] = lsqnonlin(@MO.Gdiffcompfit, parf, lb, hb, options, names, par, data(ist:ien,:));
[parf, norm] = lsqnonlin(@MO.Gdiffcompfit, parf, lb, hb, options, names, par, data);
for i = 1:length(names)
par = setfield(par, names{i}, parf(i));
end
......@@ -192,19 +272,30 @@ classdef absFcsmodel < handle
fileID = fopen(fname);
% read whole file at once
C = textscan(fileID, '%s', 'delimiter','\n');
% find Corr values
rawDataName = find(strncmp(C{1}, 'RawData =', 9)); %this reads all entries independenty of Ch1 or Ch2
% find number of channels pts and repetions
fclose(fileID);
%%
% find repetion position and channel
RPC = [];
for i=1:length(rawDataName)-1
[tok] = regexp(C{1}(rawDataName(i)),'.+_R(?<rep>(\d+))_P(?<pnt>(\d+))_K\d_Ch(?<ch>\d).raw', 'names');
tok{1}.rep = str2num( tok{1}.rep);
tok{1}.pnt = str2num( tok{1}.pnt);
tok{1}.ch = str2num( tok{1}.ch);
RPC = [RPC;struct2array(tok{1})];
positionIdx = find(strncmp(C{1}, 'Position = ', 11));
repetitionIdx = find(strncmp(C{1}, 'Repetition = ', 13));
channelIdx = find(strncmp(C{1}, 'Channel = ', 10));
for i=1:length( positionIdx)-1
tmp = C{1}(positionIdx(i));
[tokp] = regexp(tmp{:},'(\d+)', 'tokens');
tmp = C{1}(repetitionIdx(i));
[tokr] = regexp(tmp{:},'(\d+)', 'tokens');
tmp = C{1}(channelIdx(i));
[tokc] = regexp(tmp{:},'(\d+)', 'tokens');
%invert channel for compability with FA this assumes APD!!!
if str2num(char(tokc{:}))==2
ch = 1;
elseif str2num(char(tokc{:}))==1
ch = 2;
end
RPC = [RPC; str2num(char(tokr{:}))+1 str2num(char(tokp{:}))+1 ch ];
end
% index for measurements is based on repetition or measured points
idxR = 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