%%
%
% skript je převzatý z kodu viz. MÉZL, M a R JIŘÍK. Srovnání optimalizačních algoritmů pro ultrazvukovou perfúzní analýzu metodou Bolust & Burst. Ústav biomedicínského inženýrství, VUT v Brně, Ústav přístrojové techniky AV ČR, v.v.i., 5
% ve skriptu je kod pro výpočet SOMA převzatý z TVRDÍK, Josef. Evoluční algoritmy. Skripta, Přírodovědecká fakulta, Ostravská univerzita, 2004.  
clear all
clc
close all
%%
%load data_nove/sim_2logn_100times_SNR50_inp_tis_150813.mat
load('C:\škola\odevzdání\data (3)\SOMA\data_nove\sim_2logn_100times_SNR50_inp_tis_150813.mat');
% definition of settings for perfusion analysis
settings=struct(); % initiate
settings.aifmodel='2logn'; % only one model now
settings.deconvtype='blind'; % only blind is done
settings.aifinit=100; % number of initial estimates
settings.tis_scale=0; % tissue scaling - for simulated data must be 0 otherwise parameters of AIF are different than set
settings.initpar='...'; % not implemented yet
settings.disp_results='on';    % plot for each curve in tissue
% settings.display='full'; % not yet 
%[tissue,info]=us_perfan_v031pso(tissue,info,settings); % call the function
%[tissue,info]=us_perfan_v031(tissue,info,settings);
%[tissue,info]=us_perfan_v031soma(tissue,info,settings);
info.codewords{size(info.codewords,2)+1}='par';
%[filename, info]=name_file(info);    % Michal`s function
%save(filename,'data1','data2','tissue','sumimg','info');

%vysledky=cell(40,10);
% casy=cell(40,10);
% minikriteria=cell(40,10);
% kritpuvod=cell(40,10);

% function for us perfusion analysis
% IN: tissue = data, info = parameters, settings = for perfan
% OUT: tissue = result, info = updated info
% M. Mezl, 05-2015, v00 - test to make perfan in single file
% M. Mezl, 08-2015, v01 - 2gam model added for testing if concept is possible
% M. Mezl, 08-2015, v02 - in tissue now is joint vector for aif parameters  - aifpar (column 9), and TRF parameters - trfpar (column 10) - for better work with other AIF models
% M. Mezl, 09-2015, v03 - minor change - now initial estimate of AIFpar is random - from saved aifrand.mat 
%load aifrand.mat;
rois=size(tissue,1); % number of ROI in datafile
edge_start=1; % temporal = 1 for simulated data - then could be selected by user or defined automatically
flashstart=info.acq.flashstart;   % position of burst-start
flashend=info.acq.flashend;     % position of burst-end
samples=length(tissue{1,1});        % length of the tissue
flashstart=flashstart-edge_start+1; % ignoring the precontrast phase
flashend=flashend-edge_start+1;     % shift of flash indices
samples=samples-edge_start+1;       % number of samples for convolution - length of AIF
Ts=info.acq.Ts;                     % sampling period
% framelabels and work with it
if isfield(info,'framelabels')
        framelabels=info.framelabels';
else
        framelabels=zeros(size(tissue{1,1}));
end
framelabels=framelabels(edge_start:end);
tmp = ones(size(framelabels));
tmp(framelabels==-1) = 0;
framelabels = tmp;
%h=waitbar(0,'Please wait... ');

for i=1:1 %rois
    i
  
    tic;
    aifinit=settings.aifinit;
    tis=tissue{i,1}(edge_start:end);
    sum_tis=sum(tis(1:flashstart));
    if settings.tis_scale==1
        tis=tis/sum_tis;
    end
    tis_b=tis(edge_start:flashstart);
    tis_r=tis(flashend:samples);
    frm_b=framelabels(edge_start:flashstart); % XX - same for all ROIs
    frm_r=framelabels(flashend:samples); % XX - same for all ROIs
    weights=[1,1];
    N=length(tis); % XX - same for all ROIs
    Nb=length(tis_b); % XX - same for all ROIs
    NNb=2^ceil(log2((2*Nb-1))); % XX - same for all ROIs
    Nr=length(tis_r); % XX - same for all ROIs
    NNr=2^ceil(log2((2*Nr-1))); % XX - same for all ROIs
    N_cel={N,Nb,NNb,Nr,NNr};
    t_b=(1:Nb)*Ts;
    t_b_trf=t_b-Ts;
    t_r=(N-Nr+1:N)*Ts;
    t_r_trf=(1:Nr)*Ts-Ts;
    t=(1:N)*Ts;
    t_cel={t_b,t_b_trf,t_r,t_r_trf,t};
    % estimation from refill skipped 
    % now optimize
    switch settings.aifmodel
        case '3gam'
            aifsize=7;
            myfun=@(par)bab_crit_3exp_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
            parpuv0=[2.7786 .0091 .3793 .1664 .8692 .0002 .0944 .2]; % new datasets - august/september 2015 3gam
        case '2gam'
            aifsize=6;
            myfun=@(par)bab_crit_2exp_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
            parpuv0=[1.5297 0.0053 0.0618 1.9971 .1528 .5215 .2]; % new datasets - august/september 2015 2gam
        case '2logn'
            aifsize=6;
            myfun=@(par)bab_crit_2logn_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
            parpuv0=[6.9913 3.7209 0.5688 3.5039 1.9756 .7888 .2];
        case '3logn1A'
            aifsize=7;
            myfun=@(par)bab_crit_3logn1A_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
            parpuv0=[3.2495 1.9092 .7594 3.3546 .5387 3.8216 .3926 .2];
        case '3logn1sig'
            aifsize=7;
            myfun=@(par)bab_crit_3logn1sig_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
            parpuv0=[7.4203 3.7169 .5928 1.2080 1.3654 2.0085 2.2006 .2];
        otherwise 
            disp('Error is awaiting....')
    end
  %% vytvoření původní funkce 
    aif_b2=aif(1:Nb);
    aif_r2=aif(N-Nr+1:end);
    trfpar2=0.2;
    [konvol_b,konvol_r]=make_CON_1expTRF(trfpar2,aif_b2,aif_r2,t_b_trf,t_r_trf,NNb,NNr);
    model2 = [ zeros(1,edge_start-1)  konvol_b zeros(1,flashend-flashstart-1) konvol_r];
    if settings.tis_scale==1
        model2=sum_tis*model2;
    end
  %%
% SOMA algoritmus
% parametry 
popsize=1000;                   % velikost populace (smečky)          
d=aifsize+1;                    % dimenze prohledávaného prostoru 
y=tis;                          % původní­ krivka
mass=2.5;%2.5                   % přeskočení vůdce <1.1,3>
step=1.4;%1.4;                  % velikost skoku <0.11,mass>
prt=0.4;                        % pravdepodobnost vyberu dimenze <0,1>
my_eps=0.001; 
max_evals=15000;
% a=zeros(1,d);                       % vektor dolních mezí proměnných
 
a=[3 zeros(1,d-1)];
b =[7.5 4*ones(1,d-2) 1];
%b(1,2:end)=4;                   % vektor horních mezí proměnných
% xmin=a(1,1);
% xmax2=b(1,2);
% xmax=b(1,1); 

% generovani pocatecni populace
    P=zeros(popsize,d+1);
    matice=rand(popsize,d); % matice o rozmerech pocet jedincu x 1
    
    for ut=1:d
    par(:,ut)=a(1,ut)+matice(:,ut)*(b(1,ut)-a(1,ut));% nahodne generovane cislo z prislusneho rozsahu
    end
    vel=par;
    P=par;
    

       for ii=1:popsize
       nejbest=P(ii,1:d);
       aifpar=nejbest(1:aifsize);
       trfpar=nejbest(aifsize+1);
       funkce=AIFmodel22(aifpar,t,settings.aifmodel,Nb,N,Nr,trfpar,t_b_trf,t_r_trf,NNb,NNr,edge_start,flashend,flashstart,settings.tis_scale);
       %funkce=modeltyp(settings.aifmodel,nejbest,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
       krit = sum((y-funkce).^2)/(length(y));
       odhady(ii,:)=funkce;
       P(ii,d+1)=krit;
       end
       
  %vykreslení
%     figure(3)
%     plot(odhady(1,:))             %zobrazení všech odhadů první populace
%       for uz=2:popsize
%           hold on 
%           plot(odhady(uz,:))
%       end
    
  P=sortrows(P,d+1);                % seřazení jedinců podle kvality 
  pocet_kroku=fix(mass/step);       % počet skoků jednice
  func_evals=popsize;               % počet vycíslení funkce
  generace=1;                       % migrační kolo
  
  %% Algoritmus 
  tic;
while ((P(popsize,d+1)-P(1,d+1))>my_eps)&(func_evals<d*max_evals)
    leader = P(1,1:d);                  % pozice vůdce
    for il = 2:popsize                   % pro každého jednice kromě vůdce
        r0 = P(il,1:d);                  % souradnice i-teho jedince 
        delta = (leader-r0)*step;       % délka skoku jedince
        d=length(r0);                   % dimenze
        tochange=find(rand(d,1)>prt);   % vektor indexu
        
        if length(tochange)==0          % tak aslespon jednu změnu
            tochange=1+fix(d*rand(1));  % index zmeny 
        end
        
        C=zeros(pocet_kroku,d+1);       % ukladani po kazdem skoku
        r=r0;                           % souradnice jedince 
        for ii=1:pocet_kroku
            r(tochange)=r(tochange)+delta(tochange);
            r=zrcad(r,a,b);
            aifpar=r(1:aifsize);
            trfpar=r(aifsize+1);
            funkce=AIFmodel22(aifpar,t,settings.aifmodel,Nb,N,Nr,trfpar,t_b_trf,t_r_trf,NNb,NNr,edge_start,flashend,flashstart,settings.tis_scale);
            krit = sum((y-funkce).^2)/(length(y));
            %funkce=fcn(r(1,1),r(1,2),r(1,3),xx);
            fr = krit;  %modeltyp(settings.aifmodel,r,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
            C(ii,:) = [r,fr];
        end
        [fmin,indmin]=min(C(:,d+1));    % nejlepší jedinec
        yfy=C(indmin,:);                % souradnice nejlepšího 
        func_evals=func_evals+pocet_kroku;
        if yfy(d+1)<P(il,d+1)            % bod y kvalitnější 
            P(il,:)=yfy;
        end
    end % for i=2:N ... nova generace, končí migrační kolo
    generace=generace+1;
    P=sortrows(P,d+1);
end % 
%%% ulozeni vysledku
func_evals;
x_star=P(1,1:d);
fn_star=P(1,d+1)
% zobrazeni
%figure(20)
aifpar=x_star(1:aifsize);
trfpar=x_star(aifsize+1);
    switch settings.aifmodel
        case '3gam'
            aif_final=AIF_3gam(aifpar,t); % get AIF
        case '2gam'
            aif_final=AIF_2gam(aifpar,t);
        case '2logn'
            aif_final=AIF_2logn(aifpar,t);
        case '3logn1A'
            aif_final=AIF_3logn1A(aifpar,t);
        case '3logn1sig'
            aif_final=AIF_3logn1sig(aifpar,t);
        otherwise
            disp('Error is awaiting....')
    end
    aif_b=aif_final(1:Nb);
    aif_r=aif_final(N-Nr+1:end);
    [konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
    model = [ zeros(1,edge_start-1)  konvol_b zeros(1,flashend-flashstart-1) konvol_r];
        if settings.tis_scale==1
            model=sum_tis*model;
        end
    figure(4);subplot(211);plot(model(1,1:flashstart));legend('AIF model'),hold on;plot(tis(1,1:flashstart));legend('AIF');
    subplot(212);plot(model);hold on;plot(tis,'g');legend('Model','Measured');hold off;
    toc
    
%     vysledky{i,aqw}=x_star;
%     casy{i,aqw}=toc;
%     minikriteria{i,aqw}=fn_star;
%     kritf=sum((model2-model).^2)/(length(model))
%     kritpuvod{i,aqw}=kritf;
    
end
%%
%     % 2 gam
%    % tissue{i,4}=1/trfpar(1);   % mtt
%     tissue{i,5}=0;              % v - not used now
%     tissue{i,6}=0;              % f - not used now
%     tissue{i,7}=aif_final;
%     trf_temp=exp(-trfpar(1)*t_b_trf);
%     trf_temp=trf_temp/sum(trf_temp);
%     tissue{i,8}=trf_temp;
%     tissue{i,9}=aifpar;
%     tissue{i,10}=trfpar;
%     % upravit
%     %tissue{i,11}=kritbest; % rss - residual sum of squares
%     tissue{i,12}=model;
%     %tissue{i,13}=tel;
%     % now update the tissue 
%     % perfusion parameters: MTT, Fb, Vb, aifpars... (7), AIF, model,
%     % chyby...
%    % waitbar(i/rois,h);
%close(h);
info.perfan.edge_start=edge_start;
info.perfan.flashstart=flashstart;
info.perfan.flashend=flashend;
info.perfan.samples=samples;
%info.perfan.algorithm=optimget(options,'Algorithm');
info.perfan.aifmodel=settings.aifmodel;
info.perfan.aif_init=aifinit;
% 3gam
% info.tissue.descr={'Curve [-]'	'Mask'	'Mask coordinates'	'MTT [s]'	'Vb [ml/100ml]'	'Fb [ml/min/100ml]'	'AIF'	'TRF'	'beta_aif'	'A'	'tau1'	'B'	'tau2'	'C'	'tau3'	'RSS Error' 'Model' 'Elapsed time [s]'};
% 2gam
info.tissue.descr={'Curve [-]'	'Mask'	'Mask coordinates'	'MTT [s]'	'Vb [ml/100ml]'	'Fb [ml/min/100ml]'	'AIF'	'TRF'	'aifpar' 'trfpar' 'RSS Error' 'Model' 'Elapsed time [s]'};
disp('Done.')
%end

%% Funkce 

function aif=AIF_3gam(aifpar,t)
% function for computing the AIF from parameters
% IN: parameters, t
% OUT: AIF
beta_aif=aifpar(1);
A=aifpar(2);
tau1=aifpar(3);
B=aifpar(4);
tau2=aifpar(5);
C=aifpar(6);
tau3=aifpar(7);
s1=t.^beta_aif.*A.*(exp(-tau1*t));
s2=t.^beta_aif.*B.*(exp(-tau2*t));
s3=t.^beta_aif.*C.*(exp(-tau3*t));
aif=s1+s2+s3;
end
function aif=AIF_2gam(aifpar,t)
% function for computing the AIF from parameters
% IN: parameters, t
% OUT: AIF
aif=t.^aifpar(1).*aifpar(2).*exp(-aifpar(3)*t)+t.^aifpar(4).*aifpar(5).*exp(-aifpar(6)*t);
end
function aif=AIF_2logn(aifpar,t)
% function for computing the AIF from parameters
% IN: parameters, t
% OUT: AIF
aif=aifpar(1)*lognpdf(t,aifpar(2),aifpar(3))+aifpar(4)*lognpdf(t,aifpar(5),aifpar(6));
end
function aif=AIF_3logn1A(aifpar,t)
% function for computing the AIF from parameters
% IN: parameters, t
% OUT: AIF
aif=aifpar(1)*(lognpdf(t,aifpar(2),aifpar(3))+lognpdf(t,aifpar(4),aifpar(5))+lognpdf(t,aifpar(6),aifpar(7)));
end
function aif=AIF_3logn1sig(aifpar,t)
% function for computing the AIF from parameters
% IN: parameters, t
% OUT: AIF
aif=aifpar(1)*lognpdf(t,aifpar(2),aifpar(3))+aifpar(4)*lognpdf(t,aifpar(5),aifpar(3))+aifpar(6)*lognpdf(t,aifpar(7),aifpar(3));
end
function [konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr)
% function for computing convolution of AIF and TRF
% IN: parameteres, aifs, time axes, NN for FFT
% OUT: convolution result for both bolus and replenishment part
trf0=exp(-trfpar(1)*t_b_trf);
trf_b=trf0/sum(trf0);
TRF_B=fft(trf_b',NNb);
AIF_B=fft(aif_b',NNb);
KON_B=TRF_B.*AIF_B;
konvol_b=real(ifft(KON_B));
konvol_b=konvol_b(1:length(t_b_trf))';
% replenishment part
trf_r=exp(-trfpar(1)*t_r_trf);
trf_r=trf_r/sum(trf0); % CHECK this division - sum(trf0), sum(trf_r)
TRF_R=fft(trf_r',NNr);
AIF_R=fft(aif_r',NNr);
KON_R=TRF_R.*AIF_R;
konvol_r=real(ifft(KON_R));
konvol_r=konvol_r(1:length(t_r_trf))';
end
function [ssd,g]=bab_crit_3exp_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize) % HERE LETS CONTINUE""
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_3gam(aifpar,t_b);
aif_r=AIF_3gam(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [ssd,g]=bab_crit_2exp_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_2gam(aifpar,t_b);
aif_r=AIF_2gam(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [ssd,g]=bab_crit_2logn_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_2logn(aifpar,t_b);
aif_r=AIF_2logn(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [ssd,g]=bab_crit_3logn1A_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_3logn1A(aifpar,t_b);
aif_r=AIF_3logn1A(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [ssd,g]=bab_crit_3logn1sig_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_3logn1sig(aifpar,t_b);
aif_r=AIF_3logn1sig(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_2logn_blind_v00verz2(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_2logn(aifpar,t_b);
aif_r=AIF_2logn(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_2exp_blind_v00verz2(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_2gam(aifpar,t_b);
aif_r=AIF_2gam(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_3exp_blind_v00verz2(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize) % HERE LETS CONTINUE""
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_3gam(aifpar,t_b);
aif_r=AIF_3gam(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_3logn1A_blind_v00verz2(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_3logn1A(aifpar,t_b);
aif_r=AIF_3logn1A(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end
function [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_3logn1sig_blind_v00verz2(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
% criterion for optimization
%IN: parameters, tissue signals, N for fft, time axes, framelabels, weights
%OUT: now ssd, later optionally gradient
% initialization
aifpar=par(1:aifsize);
trfpar=par(aifsize+1:end);
w_b=weights(1);
w_r=weights(2);
t_b=t_cel{1};
t_b_trf=t_cel{2};
t_r=t_cel{3};
t_r_trf=t_cel{4};
NNb=N_cel{3};
NNr=N_cel{5};
% bolus and burst part of the signal
aif_b=AIF_3logn1sig(aifpar,t_b);
aif_r=AIF_3logn1sig(aifpar,t_r);
[konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
% SSD criterion 
ssd_b=sum(frm_b.*((konvol_b-tis_b).^2))*1e4;
ssd_r=sum(frm_r.*((konvol_r-tis_r).^2))*1e4;
ssd=w_b *ssd_b + w_r* ssd_r;
% gradient computing not used here...
g=1;
end

function [funkce]=modeltyp2(aifmodel,nejbest,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
    switch aifmodel
        case '3gam'
            [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_3exp_blind_v00verz2(nejbest,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '2gam'
            [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_2exp_blind_v00verz2(nejbest,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '2logn'
            [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_2logn_blind_v00verz2(nejbest,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '3logn1A'
            [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_3logn1A_blind_v00verz2(nejbest,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '3logn1sig'
            [aif_b,aif_r,konvol_b,konvol_r,ssd,g]=bab_crit_3logn1sig_blind_v00verz2(nejbest,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);

    end
    funkce = aif_b;
    
end
function [myfun]=modeltyp(aifmodel,par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize)
    switch aifmodel
        case '3gam'
            myfun= bab_crit_3exp_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '2gam'
            myfun= bab_crit_2exp_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '2logn'
            myfun= bab_crit_2logn_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '3logn1A'
            myfun= bab_crit_3logn1A_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        case '3logn1sig'
            myfun= bab_crit_3logn1sig_blind_v00(par,tis_b,tis_r,N_cel,t_cel,frm_b,frm_r,weights,aifsize);
        otherwise 
            disp('Error is awaiting....')
    end
end
%% funkce zrcad ...........................
function result = zrcad(y,a,b)
zrc=find(y<a|y>b);
for i=zrc
    while (y(i)<a(i)|y(i)>b(i))
        if y(i)>b(i)
           y(i)=2*b(i)-y(i);
        elseif y(i)<a(i)
            y(i)=2*a(i)-y(i);
        end
    end
end
result=y;
end
%% function AIF model
function model=AIFmodel22(aifpar,t,aifmodel,Nb,N,Nr,trfpar,t_b_trf,t_r_trf,NNb,NNr,edge_start,flashend,flashstart,tis_scale)
    switch aifmodel
        case '3gam'
            aif_final=AIF_3gam(aifpar,t); % get AIF
        case '2gam'
            aif_final=AIF_2gam(aifpar,t);
        case '2logn'
            aif_final=AIF_2logn(aifpar,t);
        case '3logn1A'
            aif_final=AIF_3logn1A(aifpar,t);
        case '3logn1sig'
            aif_final=AIF_3logn1sig(aifpar,t);
        otherwise
            disp('Error is awaiting....')
    end
    aif_b=aif_final(1:Nb);
    aif_r=aif_final(N-Nr+1:end);
    [konvol_b,konvol_r]=make_CON_1expTRF(trfpar,aif_b,aif_r,t_b_trf,t_r_trf,NNb,NNr);
    model = [ zeros(1,edge_start-1)  konvol_b zeros(1,flashend-flashstart-1) konvol_r];
        if tis_scale==1
            model=sum_tis*model;
        end
end
