function [aifFitted, par, SSD, residual,A_ineq,b_ineq,alb,aub,exitflag,output]=AIFfit(selected_model,t,aif0,optimAlg)

% selected model - vybrany model AIF
% tissue - tissue concentration functions - cell with tissue{k,1} being the
% k-th tissue concentration function
% aif0 - initial estimate of AIF - merena AIF
% t - time in [min]

%% Sampling period:
% Ts = t(2)-t(1); % [min]

%% Initial estimate of AIF:
[aifOptPars, alb, aub, aweights, scaleinds] = InitializeAIFParsBoundsWeights_testing_OM(selected_model);
aifNonshift = GetAIF02_OM(aifOptPars,t,selected_model);
%model scaling
scale=max(aif0)/max(aifNonshift);
aifOptParsScaled=aifOptPars;
aifOptParsScaled(scaleinds)=aifOptPars(scaleinds)*scale;

%model shifting
% aifScaled = GetAIF02_OM(aifOptParsScaled,t,selected_model);
[~,maxIdxSynt]=max(aifNonshift);
[~,maxIdxMeas]=max(aif0);
shiftSamples=maxIdxMeas-maxIdxSynt;
aifOptParsScaledShift=aifOptParsScaled;

if shiftSamples<1
    delayAIFminutes=0;
else
    delayAIFminutes=t(shiftSamples);
end
switch selected_model
    case 'Parker'
        %shifting Gauss1 Gauss2 and Sigmoid position
        aifOptParsScaledShift([3 4 10])=aifOptParsScaledShift([3 4 10])+delayAIFminutes;
    case 'ParkerShift'
        %shifting properly all signak
        aifOptParsScaledShift(end)=delayAIFminutes;
end

% fig=figure;
% ax=subplot(1,1,1,'parent',fig);
% hold(ax,'on');
% plot(ax,t,aif0,'k',t,aif,'r');
% plot(ax,t,aifNonshift,'r');
% plot(ax,t,GetAIF02_OM(aifOptParsScaledShift,t,selected_model),'g');

% optimAlg='lsqnonlin';
% optimAlg='fmincon';
% optimAlg='lsqnonlin+fmincon';


% linear inequalities
A1_A2xscale=1;
delayT1T2=0.5;
sigma2graterThenSigma1=1;
delayT2sigmoid=0.5;
expLowerThanA1=1;
switch selected_model
    case 'Parker'
        cA1=[-1 A1_A2xscale zeros(1,8)];cb1=0;%A1 is at least Xx greater comp to A2;
        cA2=[0   0 1 -1 0 0 0 0 0 0 ];cb2=-delayT1T2;%T2 je vetsi nez T1 alespon o delayT1T2
        cA3=[0   0 0  0 sigma2graterThenSigma1 -1 0 0 0 0 ];cb3=0;%sigma2 > sigma1
        cA4=[0   0 0  1 0 0 0 0 0 -1 ];cb4=-delayT2sigmoid; %sigm je posunuta doprava od T2 alespon o delayT2sigmoid
        cA5=[-expLowerThanA1 0 0 0 0 0 1 0 0 0 ];cb5=0;%amp. exponencialy je x nasobne mensi nez A1
    case 'ParkerShift'
        cA1=[-1 A1_A2xscale zeros(1,9)];cb1=0;%A1 is at least Xx greater comp to A2;
        cA2=[0   0 1 -1 0 0 0 0 0 0 0];cb2=-delayT1T2;%T2 je vetsi nez T1 alespon o delayT1T2
        cA3=[0   0 0  0 sigma2graterThenSigma1 -1 0 0 0 0 0];cb3=0;%sigma2 > sigma1
        cA4=[0   0 0  1 0 0 0 0 0 -1 0];cb4=-delayT2sigmoid; %sigm je posunuta doprava od T2 alespon o delayT2sigmoid
        cA5=[-expLowerThanA1 0 0 0 0 0 1 0 0 0 0];cb5=0;%amp. exponencialy je x nasobne mensi nez A1
end
% all
%         A_ineq=[cA1;cA2;cA3;cA4;cA5];
%         b_ineq=[cb1;cb2;cb3;cb4;cb5];
% just delays

% A_ineq=[cA1;cA2;cA3;cA4];
% b_ineq=[cb1;cb2;cb3;cb4];

% A_ineq=[cA4];
% b_ineq=[cb4];
%  delays, A1A2 amplitudes
% A_ineq=[cA1;cA2;cA4];
% b_ineq=[cb1;cb2;cb4];
A_ineq=[cA2;cA4];
b_ineq=[cb2;cb4];


%obecne omezeni kde se maji prvky vyskytnout./muze byt dale prepsano
aub([3 4 10])=t(end)/2;
switch selected_model
    case 'Parker'
        %nothing to do :)
    case 'ParkerShift'
        aub(11)=t(end);
end

switch optimAlg
    case 'doplneniG2'
        % odhad Parkera jen s jednou Gaussovkou
        par = aifOptParsScaledShift;
        par = par ./ aweights;
        alb0 = alb ./ aweights;
        aub0 = aub ./ aweights;
        
        x_start=par;
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  1  0  1  0      1      0     0    0 0]);
        x_start(x_fixed)=[0 0 Inf];
        
        optionsAIFlsqnonlin=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','trust-region-reflective','Diagnostics','off',...
            'DerivativeCheck','off','TypicalX',0.1*ones(size(x_start(1,~x_fixed))));
        
        f = @(par)CritFcnParker(par,aif0',t,selected_model, aweights,x_start,x_fixed);
        [par, SSD, residual,exitflag,output] = lsqnonlin(f,x_start(~x_fixed),alb0(~x_fixed),aub0(~x_fixed),optionsAIFlsqnonlin);
        parEnd=zeros(size(x_start));
        parEnd(1,~x_fixed)=par;
        parEnd(1,x_fixed)=x_start(x_fixed);
        par=parEnd;
        par = par .* aweights;
        
        % doplneni startovnich a omezujicich podminek pro 2. gaussovku,
        % odhad komplet parkera
        T1=par(3);
        sigma1=par(5);
        
        %lb parametry
        alb2=alb;
        alb2(4)=T1+0.6; %T2
        alb2(6)=sigma1*1; %sigma2
        %alb2(7)=par(7)/2; %
        
        %ub parametry
        aub2=aub;
        aub2(4)=T1+0.7; %T2
        aub2(6)=sigma1*15;
        
        x_start=par;
        x_start([2 4 6])=aifOptParsScaledShift([2 4 6]);
        
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  0  0  0  0      0      0     0    0 0]);
        optionsAIFlsqnonlin=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','trust-region-reflective','Diagnostics','off',...
            'DerivativeCheck','off','TypicalX',0.1*ones(size(x_start(1,~x_fixed))));
        %         par = x_start ./ aweights;
        alb0 = alb2 ./ aweights;
        aub0 = aub2 ./ aweights;
        
        f = @(par)CritFcnParker(par,aif0',t,selected_model, aweights,x_start,x_fixed);
        [par, SSD, residual,exitflag,output] = lsqnonlin(f,x_start(~x_fixed),alb0(~x_fixed),aub0(~x_fixed),optionsAIFlsqnonlin);
        parEnd=zeros(size(x_start));
        parEnd(1,~x_fixed)=par;
        parEnd(1,x_fixed)=x_start(x_fixed);
        par=parEnd;
        par = par .* aweights;
    case 'lsqnonlin2x'
        par = aifOptParsScaledShift;
        par = par ./ aweights;
        alb0 = alb ./ aweights;
        aub0 = aub ./ aweights;
        
        x_start=par;
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  1  0  1  0      1      0     0    0 0]);
        x_start(x_fixed)=[0 0 Inf];
        
        optionsAIFlsqnonlin=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','trust-region-reflective','Diagnostics','off',...
            'DerivativeCheck','off','TypicalX',0.1*ones(size(x_start(1,~x_fixed))));
        
        f = @(par)CritFcnParker(par,aif0',t,selected_model, aweights,x_start,x_fixed);
        [par, SSD, residual,exitflag,output] = lsqnonlin(f,x_start(~x_fixed),alb0(~x_fixed),aub0(~x_fixed),optionsAIFlsqnonlin);
        parEnd=zeros(size(x_start));
        parEnd(1,~x_fixed)=par;
        parEnd(1,x_fixed)=x_start(x_fixed);
        par=parEnd;
        par = par .* aweights;
        
        T1=par(3);
        sigma1=par(5);
        
        alb2=alb;
        delT1T2=0.1;%minut
        alb2(4)=T1+delT1T2;
        alb2(6)=sigma1*1.5;
%         alb2(7)=par(7)/2;
        
        aub2=aub;
        aub2(4)=T1+0.7;
        aub2(6)=sigma1*15;
        
        %%
        par = aifOptParsScaledShift;
        x_start=par;
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  0  0  0  0      0      0     0    0 0]);
        optionsAIFlsqnonlin=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','trust-region-reflective','Diagnostics','off',...
            'DerivativeCheck','off','TypicalX',0.1*ones(size(x_start(1,~x_fixed))));
        par = par ./ aweights;
        alb0 = alb ./ aweights;
        aub0 = aub ./ aweights;
        
        x_start=par;
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  0  0  0  0      0      0     0    0 0]);
        
        f = @(par)CritFcnParker(par,aif0',t,selected_model, aweights,x_start,x_fixed);
        [par, SSD, residual,exitflag,output] = lsqnonlin(f,x_start(~x_fixed),alb0(~x_fixed),aub0(~x_fixed),optionsAIFlsqnonlin);
        parEnd=zeros(size(x_start));
        parEnd(1,~x_fixed)=par;
        parEnd(1,x_fixed)=x_start(x_fixed);
        par=parEnd;
        
        % par = aifOptParsScaledShift;
        par = par ./ aweights;
        
        alb0 = alb2 ./ aweights;
        aub0 = aub2 ./ aweights;
        
        x_start=par;
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  0  0  0  0      0      0     0    0 0]);
        %         x_start(x_fixed)=[0 0 Inf];
        
        optionsAIFlsqnonlin=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','trust-region-reflective','Diagnostics','off',...
            'DerivativeCheck','off','TypicalX',0.1*ones(size(x_start(1,~x_fixed))));
        f = @(par)CritFcnParker(par,aif0',t,selected_model, aweights,x_start,x_fixed);
        [par, SSD, residual,exitflag,output] = lsqnonlin(f,x_start(~x_fixed),alb0(~x_fixed),aub0(~x_fixed),optionsAIFlsqnonlin);
        parEnd=zeros(size(x_start));
        parEnd(1,~x_fixed)=par;
        parEnd(1,x_fixed)=x_start(x_fixed);
        par=parEnd;
        par = par .* aweights;
    case 'lsqnonlin'
        par = aifOptParsScaledShift;
        x_start=par;
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  0  0  0  0      0      0     0    0 0]);
        optionsAIFlsqnonlin=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','trust-region-reflective','Diagnostics','off',...
            'DerivativeCheck','off','TypicalX',0.1*ones(size(x_start(1,~x_fixed))));
        par = par ./ aweights;
        alb0 = alb ./ aweights;
        aub0 = aub ./ aweights;
        
        x_start=par;
        %                A1 A2 T1 T2 sigma1 sigma2 alpha beta s tau
        x_fixed=logical([0  0  0  0  0      0      0     0    0 0]);
        
        f = @(par)multich_freq_dom01_altAIF_minus_y_reg(par,aif0',t,selected_model, aweights,x_start,x_fixed);
        [par, SSD, residual,exitflag,output] = lsqnonlin(f,x_start(~x_fixed),alb0(~x_fixed),aub0(~x_fixed),optionsAIFlsqnonlin);
        parEnd=zeros(size(x_start));
        parEnd(1,~x_fixed)=par;
        parEnd(1,x_fixed)=x_start(x_fixed);
        par=parEnd;
        %         figure(999)
        %         plot(t,
        
    case 'fmincon'
        optionsAIFfmincon=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','interior-point','Diagnostics','off',...
            'DerivativeCheck','off');
        f = @(par)multich_freq_dom01_altAIF_minus_y_reg_ssd(par,aif0',t,selected_model, aweights);
        [par, SSD,  exitflag,output] = fmincon(f,par,A_ineq,b_ineq,[],[],alb0,aub0,[],optionsAIFfmincon);
        %       [par, SSD, residual] = fmincon(f,par,[],[],[],[],alb0,aub0,[],optionsAIFfmincon);
        residual=[];
        
    case 'lsqnonlin+fmincon'
        optionsAIFlsqnonlin=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','trust-region-reflective','Diagnostics','off',...
            'DerivativeCheck','off','TypicalX',0.1*ones(size(x_start(1,~x_fixed))));
        optionsAIFfmincon=optimset('TolX',1e-10,'TolFun',1e-12,'Display','off','MaxFunEvals',1e10,...
            'Jacobian','off','MaxIter',100,'DiffMinChange',0.01,'DiffMaxChange',0.1,...
            'LargeScale','on','Algorithm','interior-point','Diagnostics','off',...
            'DerivativeCheck','off');
        
        f = @(par)CritFcnParker(par,aif0',t,selected_model, aweights);
        [par, SSD, residual,exitflag,output] = lsqnonlin(f,par,alb0,aub0,optionsAIFlsqnonlin);
        if par(3)>par(4)% T1>T2
            par=par([2 1 4 3 6 5 7:end]);
            warndlg('exchanging gaussian Gauss1 and Gauss2');
        end
        %                 f = @(par)multich_freq_dom01_altAIF_minus_y_reg_ssd(par,aif0',t,selected_model,ones(size(aweights)));
        %                 [par, SSD,  exitflag,output] = fmincon(f,par,[],[],[],[],alb0,aub0,[],optionsAIFfmincon);
        fmc = @(par)CritFcnParkerSSD(par,aif0',t,selected_model,ones(size(aweights)));
        [par, SSD, exitflag,output] = fmincon(fmc,par,A_ineq,b_ineq,[],[],alb0,aub0,[],optionsAIFfmincon);
        if isempty(SSD)
            error('problem');
        end
end
[aifFitted] = GetAIF02_OM(par,t,selected_model);