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

% Algo importante a mencionar: la estimación acertada de la mejor
% exponencial va a depender mucho de los parámetros iniciales que le
% ingresen a la función fit(). Averiguar los mejores valores para las
% condiciones iniciales de modo que nuestro ajuste converja a una
% buenasolución no es trivial, ni computacionalmente barato.

% Tal como se indica en el reporte Nº5 del GUIAD: "Se observa que el
% coeficiente de correlación del ajuste mejora con el número de
% exponenciales pero, al mismo tiempo, la precisión de estimación de los
% parámetros individuales empeora cuando se incluyen más de tres
% exponenciales."
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

% Se cargan los datos de infectados de los primeros 40 días de la pandemia
cov=readtable('datos_COVID.csv');

x=cov.dia;
y=cov.cantPersonasConInfeccionEnCurso;

%% Comparación de ajustes de exponenenciales sobre los datos de COVID de uruguay

% Funcion 'exp1'
[aj,gof]=fit(x,y,'exp1');

figure(1)
plot(x,y,'o',x,aj.a*exp(aj.b*x),'LineWidth',1.5)
hold on

% Funciones generadas por el usuario:
fx = fittype('a*exp(b*x)');
[aj2,gof2] = fit(x,y,fx,'StartPoint',zeros(1,2));

% Se grafican los resultados del ajuste sin puntos de inicio especificados
plot(x,aj2.a*exp(aj2.b*x),'LineWidth',1.5)

% Dado que los puntos de inicio de la simulación son aleatorios, algunas
% veces puede dar un buen ajuste y otras veces puede que el ajuste sea
% malo. Para controlar la calidad del ajuste se puede comenzar con 
% condiciones nulas e intentar iterar el ajuste, usando como punto de 
% partida los datos sugeridos por el ajuste anterior, hasta que se 
% estabilice el valor de R^2 (se pueden usar otras medidas de error).

tolerancia = 1e-10; %umbral a partir del que se detiene la simulación
r2_min = 0;
aux=1;
while (gof2.rsquare<r2_min) || ((abs(gof2.rsquare)-abs(aux))>tolerancia)
aux=gof2.rsquare;
[aj2,gof2] = fit(x,y,fx,'StartPoint',[aj2.a; aj2.b]);
end

plot(x,aj2.a*exp(aj2.b*x)-2,'LineWidth',1.5)% resto 2 para que se note que son similares al ajuste dado por exp1
hold off
ylim([0,max(y)+10])
xlabel('Día desde el 13/03/2020')
ylabel('i_d(t)')
legend('Datos medidos','Ajuste dado por exp1','Ajuste de fx sin puntos iniciales','Ajuste de fx iterado')
title('Comparación de ajustes')

%% Ajustes dados por tres exponenciales
figure(2)
plot(x,y,'o','LineWidth',1.5)
hold on

[v]=fit(x,y,'exp2');
% Funciones generadas por el usuario:
fx = fittype('a*exp(b*x)+c*exp(d*x)+e*exp(f*x)');
[aj3,gof3] = fit(x,y,fx,'StartPoint',[1;v.b;1;v.d;0;0]);

% Se grafican los resultados del ajuste con puntos de inicio nulos
plot(x,aj3.a*exp(aj3.b*x)+aj3.c*exp(aj3.d*x)+aj3.e*exp(aj3.f*x),'LineWidth',1.5)

% Dado que los puntos de inicio de la simulación son aleatorios, algunas
% veces puede dar un buen ajuste y otras veces puede que el ajuste sea
% malo. Para controlar la calidad del ajuste se puede intentar iterar el
% ajuste, usando como punto de partida los datos sugeridos por el ajuste
% anterior, hasta que se estabilice el valor de R^2 (se pueden usar otras
% medidas de error).

tolerancia = 1e-5; %umbral a partir del que se detiene la simulación
r2_min = 0.5;
aux=1;
while (gof3.rsquare<r2_min) || ((abs(gof3.rsquare)-abs(aux))>tolerancia)
aux=gof3.rsquare;
[aj3,gof3] = fit(x,y,fx,'StartPoint',[aj3.a; aj3.b; aj3.c; aj3.d; aj3.e; aj3.f]);
end

plot(x,aj3.a*exp(aj3.b*x)+aj3.c*exp(aj3.d*x)+aj3.e*exp(aj3.f*x),'LineWidth',1.5)
hold off
ylim([0,max(y)+10])
xlabel('Día desde el 13/03/2020')
ylabel('i_d(t)')
legend('Datos medidos','Ajuste de fx sin puntos iniciales','Ajuste de fx iterado')
title('Comparación de ajustes')

% Puede demorar como 1h

%% Ajustes dados por cuatro exponenciales
figure(3)
plot(x,y,'o','LineWidth',1.5)
hold on

% Funciones generadas por el usuario:
fx = fittype('a*exp(b*x)+c*exp(d*x)+e*exp(f*x)+g*exp(h*x)');
% Defino límites para los valores de cada parámetro
options = fitoptions(fx);
options.Lower = -ones(8,1)*1e10;
options.Upper = ones(8,1)*1e10;
options.StartPoint = [1 aj3.b 1 aj3.d 1 aj3.f 1 0];
[aj4,gof4] = fit(x,y,fx,options);

% Se grafican los resultados del ajuste con puntos de inicio nulos
plot(x,aj4.a*exp(aj4.b*x)+aj4.c*exp(aj4.d*x)+aj4.e*exp(aj4.f*x)+aj4.g*exp(aj4.h*x),'LineWidth',1.5)

% Dado que los puntos de inicio de la simulación son aleatorios, algunas
% veces puede dar un buen ajuste y otras veces puede que el ajuste sea
% malo. Para controlar la calidad del ajuste se puede intentar iterar el
% ajuste, usando como punto de partida los datos sugeridos por el ajuste
% anterior, hasta que se estabilice el valor de R^2 (se pueden usar otras
% medidas de error).
%%
tolerancia = 1e-10; %umbral a partir del que se detiene la simulación
r2_min = 0.7;
aux=1;
while (gof4.rsquare<r2_min) || ((abs(gof4.rsquare)-abs(aux))>tolerancia)
aux=gof4.rsquare;
[aj4,gof4] = fit(x,y,fx,'StartPoint',[aj4.a; aj4.b; aj4.c; aj4.d; aj4.e; aj4.f; aj4.g; aj4.h]);
end

plot(x,aj4.a*exp(aj4.b*x)+aj4.c*exp(aj4.d*x)+aj4.e*exp(aj4.f*x)+aj4.g*exp(aj4.h*x),'LineWidth',1.5)
hold off
ylim([0,max(y)+10])
xlabel('Día desde el 13/03/2020')
ylabel('i_d(t)')
legend('Datos medidos','Ajuste de fx sin puntos iniciales','Ajuste de fx iterado')
title('Comparación de ajustes')

% Puede demorar como 1h