clear all
clc
clf

% Input Parameters
DataFile = {'Data28062018.txt'};
TimeLabel = {'time (s)'};
SignalLabel = {'{\itS}({\itt}) = {\itE} (mV)'};
TimeUnit = {'s'};
SignalUnit = {'mV'};

t_res = 1; % Time resolution
t_inputchange = 270; % Time at which the input change is done 

cal_slope = 98.9; % Slope form pH-meter calibration in %; type 100 for optical sensors
cal_offset = 11.2;% Offset form pH-meter calibration in mV; not relevant for optical sensors

ModelType = 3;    % 1: Emperical first-order LTI model
                  % 2: Emperical hyperbolic model
                  % 3: Physical film diffusion controlled model
                  % 4: Physical ion exchange controlled model
                  % 5: Physical bulk optode diffusion controled model
k_guess = 0.02; % May be adjusted according to the model type chosen

alpha_transient = 0.99;
ParTol = 0.10;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load data
Data = importdata(DataFile{1});
t = Data(:,1) - t_inputchange; % time
y = Data(:,2); % signal output

% Calculate pH values from output signal
slope = 59.16*(cal_slope/100);
CalData = [slope,cal_offset];
pHfun = @(S,c)((7*c(1)+c(2))-S)/(c(1));
pH = pHfun(y,CalData);
activity = [pH(1) pH(end)]; % Approx. pH at steady state before and after activity change 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 3: Initial curve fit %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Starting values, LTI-model function and boundaries
a0a = 10^(activity(2)-activity(1));
drift = 0;
k = k_guess; 
S1 = mean(y(end-10:end));
td = 0;

x0 = [a0a drift k S1 slope td];

ftobj = fittype('[S1 + DRIFT*(x-td) + SLOPE*log10(A0A)].* (x <= td) + [S1 + DRIFT*(x-td) + SLOPE*log10(A0A)*(exp(-(x-td)*K))].* (x > td)');
x_up = [Inf Inf Inf Inf slope t_res];
x_low = [0 -Inf 0 -Inf slope 0];

% Initial curve fit using first-order LTI model
options = fitoptions(ftobj);
options.Robust = 'On';
options.StartPoint = x0;
options.Lower = x_low;
options.Upper = x_up;
cfobj = fit(t,y,ftobj,options);

% Update parameters, statistics
x = coeffvalues(cfobj);
a0a = x(1);
drift = x(2);
k = k_guess;
S1 = x(4);
td = 0;

S0 = S1 + log10(a0a)*slope;
S_calc = cfobj(t);
Confid95 = confint(cfobj);
Stdevs1 = (x - Confid95)/1.96;

% Cut out transient region data determined by input alpha_transient
Idx(1) = find(max(diff(y)) == diff(y))-1;
Idx(2) = find(y > S0 + alpha_transient*(S1 - S0), 1 );
IdxA = (Idx(1)-(Idx(2)-Idx(1)):Idx(2)+(Idx(2)-Idx(1)));
IdxB = [(1:Idx(1)-(Idx(2)-Idx(1))),(Idx(2)+(Idx(2)-Idx(1)):length(y))]; %Complementary to IdxA


% Noise estimation
Err = y(IdxB)-S_calc(IdxB);
sigma = std(Err);
signal = abs(log10(a0a)*slope);
SN = signal/sigma;

% Control that enough data has been collected
alpha_max = (SN - 1)/SN;
idx = find(y > S0 + alpha_max*(S1 - S0), 1 );
t_alpha_max = t(idx);
if t(end) > t_alpha_max*x(3)
    disp('Steady state has been reached within given S/N-ratio')
else
    disp('Increase tmax for data collection')
end

%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 4: Curve fitting %
%%%%%%%%%%%%%%%%%%%%%%%%%

% Starting values, model functions and boundaries
x0 = [a0a drift k S1 slope td];

if ModelType == 1
    ModelTypeLabel = {'Empirical first order LTI model (i)'};
    ftobj = fittype('[S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)].* (x <= t0) + [S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)*(exp(-(x-t0)*K))].* (x > t0)');
    x_up = [(1+ParTol)*a0a drift Inf (1+ParTol)*S1 slope t_res];
    x_low = [(1-ParTol)*a0a drift 0 (1-ParTol)*S1 slope 0];

elseif ModelType == 2 
    ModelTypeLabel = {'Empirical hyperbolic model (ii)'};
    ftobj = fittype('[S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)].* (x <= t0) + [S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)*(1-(((x-t0)*K)./(1+((x-t0)*K))))].* (x > t0)');
    x_up = [(1+ParTol)*a0a drift Inf (1+ParTol)*S1 slope t_res];
    x_low = [(1-ParTol)*a0a drift 0 (1-ParTol)*S1 slope 0];
    
elseif ModelType == 3
    ModelTypeLabel = {'Physical film diffusion controlled model (iii)'};
    ftobj = fittype('[S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)].* (x <= t0) + [S1 + DRIFT*(x-t0) + SLOPE*log10(1-(1-A0A)*exp(-(x-t0)*K))].* (x > t0)');
    x_up = [(1+ParTol)*a0a drift Inf (1+ParTol)*S1 slope t_res];
    x_low = [(1-ParTol)*a0a drift 0 (1-ParTol)*S1 slope 0];
    
elseif ModelType == 4
    ModelTypeLabel = {'Physical ion exchange controlled model (iv)'};
    ftobj = fittype('[S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)].* (x <= t0) + [S1 + DRIFT*(x-t0) + SLOPE*log10(1-(1-A0A)*(1./(sqrt((x-t0)*K)+1)))].* (x > t0)');
    x_up = [(1+ParTol)*a0a drift Inf (1+ParTol)*S1 slope t_res];
    x_low = [(1-ParTol)*a0a drift 0 (1-ParTol)*S1 slope 0];
    
elseif ModelType == 5
    ModelTypeLabel = {'Physical bulk optode diffusion controlled model (v)'};
    ftobj = fittype('[S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)].* (x <= t0) + [S1 + DRIFT*(x-t0) + SLOPE*log10(A0A)*(8/pi^2)*((1/1).*exp(-1*K*(x-t0)) + (1/9).*exp(-9*K*(x-t0)) + (1/25).*exp(-25*K*(x-t0)) +(1/49).*exp(-49*K*(x-t0)))].* (x > t0)');
    x_up = [(1+ParTol)*a0a drift 0.65 (1+ParTol)*S1 slope t_res];
    x_low = [(1-ParTol)*a0a drift 0 (1-ParTol)*S1 slope 0];
end

% Curve fitting
options = fitoptions(ftobj);
options.Robust = 'On';
options.StartPoint = x0;
options.Lower = x_low;
options.Upper = x_up;
[cfobj,gof] = fit(t(IdxA),y(IdxA),ftobj,options);

% Update parameters and fit statistics
x = coeffvalues(cfobj);
a0a = x(1);
k = x(3);
S1 = x(4);
td = x(6);

S0 = S1 + log10(a0a)*slope;
S_calc = cfobj(t);
if ModelType == 5
    S_calc(isnan(S_calc)) = S1 + drift*(t(isnan(S_calc))-td) + slope*log10(a0a);
end

Confid95 = confint(cfobj);
Stdevs2 = (x - Confid95)/1.96;
SSE = sum((y - S_calc).^2);
SST = sum((y - mean(y)).^2);
Rsquare = 1 -(SSE/SST);

% Noise estimation
Err = y(IdxB)-S_calc(IdxB);
sigma = std(Err);
signal = abs(log10(a0a)*slope);
SN = signal/sigma;

% Control that robustly determined drift
dSdt = smooth(diff(y),0.10,'lowess')./diff(t);
if dSdt(end) > Stdevs1(1,2)
    disp('Drift determination is not conflicting curve fit')
else
    disp('Drift is not robustly determined')
end

% Parameters for plot
xPos = 0.25*t(end);
yPos = min(y);
yPosDelta = abs(max(y)-min(y));

t_plot = (min(t):0.01*(t(2)-t(1)):max(t));
S_plot = cfobj(t_plot);
S_plot(isnan(S_plot)) = S1 + drift*(t_plot(isnan(S_plot))-td) + slope*log10(a0a);
IdxC(1) = find(t(IdxA(1)) < t_plot, 1 );
IdxC(2) = find(t_plot < t(IdxA(end)), 1, 'last' );
IdxC = (IdxC(1):1:IdxC(2));

% Calculate t_alpha
syms w
LHside90 = 0.90*(S1-S0)+S0;
LHside95 = 0.95*(S1-S0)+S0;
LHside99 = 0.99*(S1-S0)+S0;

if ModelType == 1
    RHside = S1 + slope*log10(a0a)*exp(-(w-td)*k);
elseif ModelType == 2
    RHside = S1 + slope*log10(a0a)*(1-(((w-td)*k)./(1+((w-td)*k))));
elseif ModelType == 3
    RHside = S1 + slope*log10(1-(1-a0a)*exp(-(w-td)*k));
elseif ModelType == 4
    RHside = S1 + slope*log10(1-(1-a0a)*(1./(sqrt((w-td)*k)+1)));
end

if ModelType <= 4
    t90 = eval(solve(LHside90 == RHside,w));
    t95 = eval(solve(LHside95 == RHside,w));
    t99 = eval(solve(LHside99 == RHside,w));
elseif ModelType == 5
    t90 = t_plot(round(mean(find(abs(S_plot - LHside90) < 0.05))));
    t95 = t_plot(round(mean(find(abs(S_plot - LHside95) < 0.05))));
    t99 = t_plot(round(mean(find(abs(S_plot - LHside99) < 0.05))));
end

%%%%%%%%%%
% Figure %
%%%%%%%%%%
f1 = figure(1);
p = uipanel('Parent',f1,'BorderType','none');
p.Title = ModelTypeLabel;
p.TitlePosition = 'centertop';
p.FontSize = 14;
p.FontWeight = 'bold';
set(gcf,'color','w');

subplot(1,2,1,'Parent',p)
hold on
plot(t,y,'o','MarkerSize',3,'Color',[80 180 200]/255)
plot(t_plot,S_plot,'-','LineWidth',1,'Color',[0 0 0]/255)
plot([t(IdxA(1)) t(IdxA(1))],[yPos-0.1*yPosDelta yPos+1.4*yPosDelta],'--k')
plot([t(IdxA(end)) t(IdxA(end))],[yPos-0.1*yPosDelta yPos+1.4*yPosDelta],'--k')

text(xPos,yPos + 0.55*yPosDelta,['{\ita}_{0}/{\ita}_{inf} = ',num2str(x(1),'%1.5f'),...
    ' (',num2str(Stdevs1(1,1),'%1.5f'),')'])
text(xPos,yPos + 0.45*yPosDelta,['{\itS}_{inf} = ',num2str(x(4),'%1.2f'),...
    ' (',num2str(Stdevs1(1,4),'%1.2f'),') ',SignalUnit{1}])
text(xPos,yPos + 0.35*yPosDelta,['{\itx}_{drift} = ',num2str(x(2),'%1.1e'),...
    ' (',num2str(Stdevs1(1,2),'%1.1e'),') ',SignalUnit{1},'/',TimeUnit{1}])
text(xPos,yPos + 0.25*yPosDelta,['\sigma = ',num2str(sigma,'%1.2f'),' ',SignalUnit{1}])
title('Overall fit')
xlabel(TimeLabel)
ylabel(SignalLabel)
legend('Observed response','Curve fit','Transient region')
ylim([yPos-0.1*yPosDelta yPos+1.4*yPosDelta])
grid minor

subplot(1,2,2,'Parent',p)
xPos = 0.3*(t(IdxA(end))-x(6));
hold on
plot(t(IdxA),y(IdxA),'o','MarkerSize',3,'Color',[80 180 200]/255)
plot(t_plot(IdxC),S_plot(IdxC),'-','LineWidth',1,'Color',[0 0 0]/255)
text(xPos,yPos + 0.55*yPosDelta,['{\itk} = ',num2str(x(3),'%1.3f'),...
    ' (',num2str(Stdevs2(1,3),'%1.3f'),') ',TimeUnit{1},'^{-1}'])
text(xPos,yPos + 0.45*yPosDelta,['{\itt}_{d} = ',num2str(x(6),'%1.2f'),...
    ' (',num2str(Stdevs1(1,6),'%1.2f'),') ',TimeUnit{1}])
text(xPos,yPos + 0.25*yPosDelta,['R^2 = ',num2str(Rsquare,'%1.4f')])
text(xPos,yPos + 0.15*yPosDelta,['SSE = ',num2str(SSE,'%4.0f')])
title('Fit in transient region')
xlabel(TimeLabel)
ylabel(SignalLabel)
legend('Observed response','Curve fit')
ylim([yPos-0.1*yPosDelta yPos+1.4*yPosDelta])
grid minor