% % % % 
% funkce slouží k výpočtu spojitého GA, RMSE a Rkvadrátu
 clear all
 clc
%  
load('funkce2.mat');

fcn=@(a,b,x) 1000000*gampdf(x,a,b); 
%fcn=@(a,b,x) 1000000*lognpdf(x,a,b)
%fcn=@(a,b,x) 100000*exp(-(log(x)-a).^(2)./(2.*b.^(2)))/(x.*b.*sqrt(2*pi))
%fcn=@(a,b,x) 400000*(exp(a)./b).*(sqrt(b./(x))).*(a./(2*pi)).^(3/2).*exp((-1/2).*a.*((b./x)+(x./b)))
% fcn=@(a,b,x) 1000000*(exp(a)./b).*(sqrt((b.*a)./(x*2*pi))).*exp((-1/2).*a.*((b./x)+(x./b)))

%ukol 1
x = 0:0.01:60;
x(1,1)=0.01;
x=x(1,1:6001);
 y=funkce2;
 y=downsample(y,10);
 y=y';

% % %% ukol 5 - GA se spojitym vyjadrenim hodnot
jedincu=20; % pocet jedincu v populaci
elite=1;    % volba elitarstvi
pst_m=0.2;  % pravdepodobnost mutace
max_iter=50;   % maximalni pocet iteraci
% 
% % generovani pocatecni populace
matice=rand(jedincu,2); % matice o rozmerech pocet jedincu x 3 (mame 3 promenne)
populace(:,1)=0+matice(:,1)*(10-(0)); % nahodne generovane cislo z prislusneho rozsahu
populace(:,2)=0+matice(:,2)*(10-(0));
%populace(:,3)=-10+matice(:,3)*(10-(-10));

for ind=1:jedincu
    funkce = fcn(populace(ind,1),populace(ind,2),x); % stanoveni kvality jako krit. funkce * (-1)
    kvalita(ind,1) = -sum((y-funkce).^2)/(length(y)); %hledame min. -> zaporna hodnota
end
% plot(populace(:,1),populace(:,2),'rx'); % vykresleni pocatecni populace
stagnace=0;
fbest=-Inf;
ite=1;
%%
%% samotny algoritmus
while (ite<=max_iter)&& (stagnace<30) % stejne jako minule - podm. max poctu iteraci a stagnace reseni
    for ind=1:jedincu
        funkce = fcn(populace(ind,1),populace(ind,2),x); % stanoveni kvality jako krit. funkce * (-1)
        kvalita(ind,1) = -sum((y-funkce).^2)/(length(y)); %% stanoveni kvality jako krit. funkce * (-1)
    end
    [m, pos]=max(kvalita);
    elit=populace(pos,:);
    % selekce a reprodukce - vyberu lepsi polovinu jedincu a ty se budou
    % reprodukovat - vystupem je populace se stejnym poctem clenu
    populace=repr(populace,kvalita,fcn,y);
    % mutace - nekteri jedinci jeste podlehnou mutaci
    populace=mutace_s(populace,pst_m,ite/max_iter);
    if elite==1
        populace(1,:)=elit;
    end
    for ind=1:jedincu
        funkce = fcn(populace(ind,1),populace(ind,2),x); % stanoveni kvality jako krit. funkce * (-1)
        kvalita(ind,1) = -sum((y-funkce).^2)/(length(y)); % stanoveni kvality jako krit. funkce * (-1)
    end
    if abs(max(kvalita)-fbest)<1e-6 % pokud se bude nejlepsi jedinec lisit o 1e-6 od nejlepsiho v minulem kroku
        stagnace=stagnace+1;
    else
        fbest=max(kvalita);
        stagnace=0;
    end
    pause(0.1)
    prubeh_kvality(ite)=mean(kvalita);
    prubeh_elity(ite)=max(kvalita);
    [bes pos] = max(kvalita);
    disp(['Nejlepsi jedinec: [' num2str(populace(pos,1)) ', ' num2str(populace(pos,2)) ']. Stagnace: ' num2str(stagnace)]);
    ite=ite+1;
end

figure(4)
plot(prubeh_kvality); hold on;
plot(prubeh_elity,'r');
title('Prubeh hledani extremu');
xlabel('Iterace');
ylabel('Kvalita');
legend('Prumerna kvalita','Elitni jedinec','Location','southeast');
disp('______________________________________________');
disp(['Reseni: [' num2str(populace(pos,1)) ', ' num2str(populace(pos,2)) '].']);

vysl2 = fcn(populace(pos,1),populace(pos,2),x);
figure(5)
plot(x,y,'Color','b');
hold on;
plot(x,vysl2,'Color','r');
% % % výpočet RMSE a R^2
 RMSE= sqrt(mean((y - vysl2).^2)) 

 
 % Root Mean Squared Error
R = corrcoef (y,vysl2);
Rsq = R (1,2).^2