% liftinglockdown.m
% Script to evaluate what could happen once lockdown is lifted based on different R0 values
%   by default tau is assumed to be equal to 10 days
%   current script explore the consequence around Paris
%   RESULTS ARE PROVIDED "AS IS"
%
%   INRAe\Olivier Vitrac - 2020-01-05

% LOAD RAW DATA: selection of departments arount Paris
outputfolder = 'after_lifting_lockdown/';
if ~exist(outputfolder,'dir'), mkdir(outputfolder), end
if ~exist('R','var') || isempty(R)
    R=COVID19Modelingv2('ParisSuburb'); % the current study focuses around Paris
end
% sort entities
nR = length(R);
[~,ind] = sort([R.entity],'ascend');
R = R(ind);
% add population
% Source 2017: https://fr.wikipedia.org/wiki/Liste_des_d%C3%A9partements_fran%C3%A7ais_class%C3%A9s_par_population_et_superficie
% Source 2020: https://www.ined.fr/fr/tout-savoir-population/chiffres/france/structure-population/regions-departements/
pop = {
    "DPT 75,Paris"               2187526     
    "DPT 77,Seine-et-Marne"      1403997
    "DPT 91,Essonne"             1296130
    "DPT 92,Hauts-de-Seine"      1609306
    "DPT 93,Seine-Saint-Denis"   1623111
    "DPT 94,Val-de-Marne"        1387926
    "France"                     64897954
    "Région Île-de-France"       12278210
};
for i=1:nR, R(i).population = pop{[pop{:,1}]==R(i).entity,2}; end
% additional definitions
datefmt = 'dd/mm/yyyy';
E = datenum('11/05/2020',datefmt);  % end of lockdown date
F = E + 100;                         % end of predictions
nEF = F - E;
col = [tooclear(cbrewer('qual','Set2',nR-2)); 0 0 0; .5 .5 .5];
lw = ones(nR,1) * 2; lw(end-1:end) = lw(end-1:end)+2;
project = ['Paris_lift_' datestr(today,'yyyy-mm-dd')];
removelast = @(x) [x(1:end-1);{''}];
removelastyticklabel = @(h) set(h,'yticklabel',removelast(get(h,'yticklabel')));
close all

%% SCENARIO after the lockdown is lifted: date E
% simple model: n = n0 * R ^ (t/tau)
% Source: https://en.wikipedia.org/wiki/Basic_reproduction_number
% Reprudction number: R0
%R0 = 1.0;  % Reproduction number (higher than the German value, assuming worst case)
tau = 10; % days  % contamination period, no latence

for R0 = [1:.01:1.05 1.1:.05:1.3]

    for i=1:nR
        j = find( R(i).data.date == E );
        newinfected = zeros(1,nEF);
        buffer = R(i).data.pdfpred(j+(-tau+1:0));
        for k = 1:nEF
            newinfected(k) = sum(buffer*R0/tau);
            buffer(1:end-1) = buffer(2:end);
            buffer(end) = newinfected(k);
        end
        R(i).data.iE = j;
        R(i).data.tnewinfected = R(i).data.date(j)+(0:nEF);
        R(i).data.newinfected = [R(i).data.pdfpred(j) newinfected];
    end
    
    
    %% EXTRAPOLATION BASED ON CURRENT EVOLUTION (no change)
    clf
    formatfig(gcf,'figname',sprintf('%s_R0=%0.2f',project,R0),'position',[956 312 730 950],'paperposition',[1.5000    4.1750   18.0000   21.3500])
    hs = subplots(1,[.4 1],0,0); hax=axes('position',[.4936 .2846 .3813 .20]);
    hp = zeros(nR,1);
    tickfmt = 'mmm,dd';
    for i=1:nR
        population = R(i).population;
        j = R(i).data.iE;
        newcases = R(i).data.ypredCI(j,2) + [0 cumsum(R(i).data.newinfected(2:end))];
        
        subplot(hs(1)), hold on
        plot(R(i).data.date(1:j),1000*R(i).data.ypredCI(1:j,2)/population,'-','color',col(i,:),'linewidth',4)
        plot(R(i).data.tnewinfected,1000*newcases/population,'-','color',col(i,:),'linewidth',2)
        plot(R(i).data.sampledDate,1000*R(i).data.infected/population,'o','markerfacecolor',col(i,:),'markeredgecolor','w','markersize',5)
        
        for cax = [hax hs(2)]
            subplot(cax), hold on
            hp(i) = plot(R(i).data.date(1:j),1000*R(i).data.pdfpred(1:j)/population,'-','color',col(i,:),'linewidth',lw(i));
            ind = find((R(i).data.date>=today) & (R(i).data.date<=E));
            if (cax==hax) && (R(i).entity=="DPT 75,Paris")
                text(R(i).data.date(ind(1)),1000*R(i).data.pdfpred(ind(1))/population,'\leftarrow TODAY','fontsize',12,'horizontalalignment','left','verticalalignment','middle','rotation',270,'fontweight','bold')
                text(R(i).data.date(ind(end)),1000*R(i).data.pdfpred(ind(end))/population,'\leftarrow May, 11^{th}','fontsize',12,'horizontalalignment','left','verticalalignment','middle','rotation',90,'fontweight','bold')
            end
            plot(R(i).data.date(ind),1000*R(i).data.pdfpred(ind)/population,'-','color',col(i,:),'linewidth',2/3*lw(i))
            plot(R(i).data.tnewinfected,1000*R(i).data.newinfected/population,'-','color',col(i,:),'linewidth',lw(i)/2)
        end
    end
    title(hs(1),sprintf('Lockdown lifted on %s followed by \\fontsize{12}R_0\\approx\\bf%0.3g\\rm,\\tau=%d days',datestr(E,'yyyy-mm-dd'),R0,tau),'fontsize',12,'fontweight','normal')
    ylabel(hs(1),'Total cases per 1000','fontsize',14)
    ylabel(hs(2),'Daily cases per 1000','fontsize',14)
    ylabel(hax,'Daily cases per 1000','fontsize',10)
    
    
    xlabel(hs(2),'date','fontsize',14)
    formatax(hax,'fontsize',10)
    datetick(hs(1),'x',tickfmt)
    datetick(hs(2),'x',tickfmt)
    datetick(hax,'x',tickfmt)
    formatax(hs,'fontsize',14,'xlim',[datenum('20/03/2020',datefmt) datenum('20/07/2020',datefmt)],'yscale','linear')
    set(hax,'xlim',[datenum('20/04/2020',datefmt) datenum('11/06/2020',datefmt)])
    subplot(hs(2))
    legend(hp,regexprep([R.entity],'^DPT ',''),'location','best','fontsize',10,'box','off')
    removelastyticklabel(hs(2))
    
    % print
    print_png(300,[get(gcf,'filename') '.png'],outputfolder,'',0,0,0)
    print_pdf(300,[get(gcf,'filename') '.pdf'],outputfolder,'nocheck')
end % next R0