% qualitative multiscale corona random walk model
% author: stefan luding, 01-April-2020
%
clear all
close all
format compact

warning('on','MATLAB:RandStream:ActivatingLegacyGenerators')
warning('on','MATLAB:RandStream:ReadingInactiveLegacyGeneratorStates')
sd=101; % rng seed
srng0 = rng(sd) % random number settings

% which type of run 4 different runs => 4 plots
irun=0; % compare effect of different infection rates
irun=1; % compare effect of different flying rates
irun=2; % compare effect of different random step sizes
irun=3; % compare effect of different cell-sizes
irun>10; % compare effect of duration of different actions (see below)
% set one irun value, or run through all
auto_mode=0; % use this to ask for printing 
%for irun=[0 1 2 3 11 12 13 14 16 17]
%for irun=[11 12 14 16 17]
for irun=[27 12 11]
    auto_mode=1 % do not ask - print automatically
    irun
    figure(10); clf;
    figure(11); clf;
    figure(12); clf;

% time definitions
dt=0.2; % every calculation loop
itime=500; %00 % units of days/dt
ipstep=itime/20; % output every ipstep time steps

dayaction1=200 % day when action is taken
dayaction2=200 % day when action is stopped

action1=0;
if(irun > 10)
    action1=irun-10 % actions can be set together binary mode - e.g. action=7 activates 1,2,4
end
if(irun > 20)
    action1=irun-20 % actions can be set together binary mode - e.g. action=7 activates 1,2,4
end
if(bitand(1,action1)>0) % just stop practically all long-distance flights by 0.01
    ffly=0.01
end
if(bitand(2,action1)>0) % lock-down, immobilize 0.9 of the population
    fslow=0.1
end
if(bitand(4,action1)>0) % reduce infection rate by factor 0.5
    fifrate=0.5
end

%% space definitions
Dim=2; % dimension 
Lsys=100;
Lsys=2000 % system size Lsys x Lsys
% setup fields
%N=300; % population
N=500000 % population
%
%p0=0.5e-5; % initial probability for infection
p0=2e-5; % initial probability for infection
%
vifrate=[.001 .002 .004 .008] *dt; % infection rate/probability
%vifrate=[.002 .004 .008 .016] *dt; % infection rate/probability
mkr=['+' 'x' 'o' 'square' 'diamond'];
%vifrate=[.001]; % infection rate/probability
sdays = 20; % days of sickness
stx = 0.6; % incubation period => after 1-stx of the sickness infectious
rhealing=dt/sdays; % healing rate = 1 / days of infection
mort=0.1; % mortality fraction
rdeath=mort*dt/sdays; % tested every timestep during days of sickness
%
idist=20 % critical distance
imax=floor(Lsys/idist) % cell-size within range
Ncell=N/(imax*imax)
step=100*dt % stepsize
probfly=0.001*dt; % probability of flights per day
%stepfly=Lsys/5 % distance multiplier for flights (old)
stepfly=Lsys/2 % distance multiplier for flights (new)
%
for ifr=1:4 % loop over different settings
    % set RNG everytime same seed from above
    rng(srng0)
    rand(5,1)
    % default irun=0
    ifrate=vifrate(ifr)
    if(irun==1)
        ifrate=vifrate(3) % default
        if(ifr==1) 
            probfly=0.00025*dt %
        elseif(ifr==2)
            probfly=0.0005*dt %
        elseif(ifr==3)
            probfly=0.001*dt % default
        elseif(ifr==4)
            probfly=0.002*dt %
        end
    elseif(irun==2)
        ifrate=vifrate(3) % default
        if(ifr==1)
            step=25*dt
        elseif(ifr==2)
            step=50*dt
        elseif(ifr==3)
            step=100*dt % default
        elseif(ifr==4)
            step=200*dt
        end
    elseif(irun==3)
        ifrate=vifrate(3) % default
        if(ifr==1)
            idist=10 % critical distance
        elseif(ifr==2)
            idist=15
        elseif(ifr==3)
            idist=20 % default
        elseif(ifr==4)
            idist=25
        end
        imax=floor(Lsys/idist) % cell-size within range
    elseif(irun==4)
        ifrate=vifrate(3) % default
        if(ifr==1)
            step=25*dt
        elseif(ifr==2)
            step=50*dt
        elseif(ifr==3)
            step=100*dt % default
        elseif(ifr==4)
            step=200*dt
        end
    elseif((irun>10)&&(irun<20))
        ifrate=vifrate(3) % default
        if(ifr==1)
            dayaction1=30
            dayaction2=30
        elseif(ifr==2)
            dayaction1=30
            dayaction2=45
        elseif(ifr==3)
            dayaction1=30
            dayaction2=60
        elseif(ifr==4)
            dayaction1=30
            dayaction2=75
        end
    elseif((irun>20)&&(irun<30))
        ifrate=vifrate(3) % default
        if(ifr==1)
            dayaction1=10
            dayaction2=50
        elseif(ifr==2)
            dayaction1=20
            dayaction2=65
        elseif(ifr==3)
            dayaction1=30
            dayaction2=75
        elseif(ifr==4)
            dayaction1=40
            dayaction2=85
        end
    end
    
    % set initial state
    x00=zeros(N,Dim); % particle positions
    x00=rand(N,2)*Lsys; % initial postions
    x0=x00;
    % initial infections
    iif=(rand(N,1)<p0);
    state=zeros(N,1); % status of particle
    state(iif)=1;
    
    % reset time series for output
    infected=zeros(itime+1,1); % time-series of infected
    infected(1)=sum(state>0)/N; % initial value
    infex=zeros(itime+1,1); % time-series of infectious
    infex(1)=sum((state>0) & (state<=stx))/N; % initial value
    immune=zeros(itime+1,1); % time-series of immune
    immune(1)=0;
    ideath=zeros(itime+1,1); % time-series of dead
    ideath(1)=0;
    i0=zeros(itime+1,1); % time-series of non-infected
    i0(1)=1-infected(1)-immune(1)-ideath(1);   
    % time-series for output
    timex=zeros(itime+1,1);
    timex(1)=0;

% loop over time-steps itime
for it=1:itime
    %it
    x=x0;
    % take actions at a certain day1
    if(timex(it)+dt/2>=dayaction1 && timex(it)-dt/2<dayaction1)
        action1
        timex(it)
        if(bitand(1,action1)>0) % just stop practicall all long-distance flights by 0.01
            stepfly=stepfly*ffly
        end
        if(bitand(2,action1)>0) % lock-down, immobilize 0.9 of the population
            step=step*fslow
        end
        if(bitand(4,action1)>0) % reduce infection rate by factor 0.5
            ifrate=ifrate*fifrate
        end
    end
    % undo the action1 on day2
    if(timex(it)+dt/2>=dayaction2 && timex(it)-dt/2<dayaction2)
        action1
        timex(it)
        if(bitand(1,action1)>0) % just stop practicall all long-distance flights by 0.01
            stepfly=stepfly/ffly
        end
        if(bitand(2,action1)>0) % lock-down, immobilize 0.9 of the population
            step=step/fslow
        end
        if(bitand(4,action1)>0) % reduce infection rate by factor 0.5
            ifrate=ifrate/fifrate
        end
    end
    % identify state of p.
    istate=state>0; % all infected
    istx=state>0&state<stx; % fully infectious
    iclear=state==-1; % immune after infection
    imort=state==-2; % dead 
    izero=state==0; % non infected (yet)
    ifly=rand(N,1)<probfly; % random flyers => frequent flyers would remain same
    
    rrr=rand(N,Dim)-0.5; % get random numbers for movement
    x(imort,1:2)=x(imort,1:2)+(rrr(imort,1:2))*0; % death no move
    x(istx,1:2)=x(istx,1:2)+(rrr(istx,1:2))*step*0.2; % infected with symptoms, after incubation move 1/5
    x(izero,1:2)=x(izero,1:2)+(rrr(izero,1:2))*step; % random step
    x(iclear,1:2)=x(iclear,1:2)+(rrr(iclear,1:2))*step; % random step
    rrr=rand(N,Dim)-0.5; % get random numbers for flights
    x(ifly,1:2)=x(ifly,1:2)+(rrr(ifly,1:2))*stepfly; % random step*factor fly
    % here one could keep the same people as flying again and again   
    % periodic boundaries
    ff=(x>Lsys); x(ff)=x(ff)-Lsys;
    ff=(x<0); x(ff)=x(ff)+Lsys;
    
    if( mod(it,ipstep)==0 ) % prepare figure for population plots
        timex(it)+dt
        Ncell_infected=sum(infected)/(imax*imax);
        figure(11)
        subplot(2,2,ifr)
        %plot(x(:,1),x(:,2),'b.') % plot all non-infected 
        axis([-0.1*Lsys Lsys*1.1 -0.1*Lsys Lsys*1.1]);
    end
    % check distances (same cell)
    icell=floor(x/idist)+1;
    
    ic=zeros(imax+1,imax+1); % cells reset 
    for i=1:N
%         i
%         state(i)
        if(state(i)>0) % count to number of infectious per cell ic
            state(i)=state(i)-rhealing; % healing in idays
        end
        ic0=ic(icell(i,1),icell(i,2));
        if(state(i)>0 && state(i)<stx) % count to number of infectious per cell ic
            if(rand()<rdeath)
                state(i)=-2;
            else
                ic(icell(i,1),icell(i,2))=ic0+1;
            end
        end
    end
    %state(1:6)
    %ic
    for i=1:N
        %i
        ic0=ic(icell(i,1),icell(i,2)); % number of infected per cell
        % WARNING - this is not proper counting!!! but linear in N
        if(ic0>=1)
            %ic0;
            if(state(i)==0) 
                risk=rand();
                if(risk<ifrate*ic0)
                    state(i)=1; % infection
                end
            end
        end
        if(state(i)<0 && state(i)>-1) % healed + immunized
            state(i)=-1;
        end
    end
    %state counting -> time-series
    timex(it+1)=it*dt;
    istate=state>0;
    istx=state>0&state<stx; % fully infectious
    iclear=state==-1;
    imort=state==-2;
    izero=state==0;
    infected(it+1)=sum(istate)/N;
    infex(it+1)=sum(istx)/N;
    immune(it+1)=sum(iclear)/N;
    ideath(it+1)=sum(imort)/N;
    i0(it+1)=sum(izero)/N;
    
    % plot situation
    if( mod(it,ipstep)==0 )
        %figure(2); surf(ic)
        figure(11)
        subplot(2,2,ifr)
        hold on; 
        plot(x(istate,1),x(istate,2),'ro','MarkerSize',0.4)
        plot(x(iclear,1),x(iclear,2),'gx','MarkerSize',0.3)
        plot(x(imort,1),x(imort,2),'kv','MarkerSize',0.8)
        cpara = [ ' run=' num2str(irun) ]; %' ifrate=' num2str(ifrate/dt) ];
        title(cpara)
        axis equal
        axis([-Lsys*0.02 Lsys*1.02 -Lsys*0.02 Lsys*1.02]);
    end
    % set for next step
    x0=x;
    
    % plotting - not every time-step, but every ipsteps
    if( mod(it,ipstep)==0 )
        figure(10)
        
        nzoom=100; % zoom in-factor

        pp=1:it;
        % plot ratios linear
        subplot(2,2,1)
        semilogy(timex(pp),infected(pp),'r')
        hold on
        semilogy(timex(pp),immune(pp),'g')
        semilogy(timex(pp),ideath(pp),'k')
        semilogy(timex(pp),i0(pp),'b')
        semilogy(timex(pp),infex(pp),'m')
        semilogy(timex(pp(end)),infected(pp(end)),'r','Marker',mkr(ifr),'MarkerSize',4.5)
        %axis([0 timex(it) p0*0.8 1])
        subplot(2,2,2)
        hold on
        plot(timex(pp),infected(pp),'r')
        plot(timex(pp(end)),infected(pp(end)),'r','Marker',mkr(ifr),'MarkerSize',4.5)
        plot(timex(pp),immune(pp),'g')
        plot(timex(pp),ideath(pp),'k')
        plot(timex(pp),i0(pp),'b')
        %axis([0 timex(it) 0 1/nzoom])   
        % plot 10*zoom
        subplot(2,2,3)
        hold on
        plot(timex(pp),infected(pp)*N,'r')
        plot(timex(pp(end)),infected(pp(end))*N,'r','Marker',mkr(ifr),'MarkerSize',4.5)
        plot(timex(pp),immune(pp)*N,'g')
        plot(timex(pp),ideath(pp)*N,'k')
        plot(timex(pp),i0(pp)*N,'b')
        axis([0 timex(it) 0 N/nzoom/10])
        % plot numbers
        subplot(2,2,4)
        hold on
        plot(timex(pp),infected(pp)*N,'r')
        plot(timex(pp(end)),infected(pp(end))*N,'r','Marker',mkr(ifr),'MarkerSize',4.5)
        plot(timex(pp),immune(pp)*N,'g')
        plot(timex(pp),ideath(pp)*N,'k')
        plot(timex(pp),i0(pp)*N,'b')
        axis([0 timex(it) 0 N/nzoom])
    end
    pause(0.001)

end
    figure(10)
    subplot(2,2,1)
    if(ifr==1)
        cl1 = [string('infected'); string('immune'); string('dead'); string('no infection'); string('exp(t/10)') ];
        cl = [string('infected'); string('immune'); string('dead'); string('no infection') ];
        plot(timex,p0*exp(timex/10),'y--','LineWidth',2)
    end

    xlabel('days')
    ylabel('fraction')
    title('log-plot of fractions')
    axis([0 timex(end) 1e-6 1])
    yticks([1e-5 1e-4 1e-3 1e-2 0.1 1])
    if(ifr==1)
        legend(cl1,'Location','northwest')
    end
    %
    subplot(2,2,2)
    xlabel('days')
    ylabel('fraction')
    title('linear-plot of fractions')
    axis([0 timex(end) 0 1/nzoom])
    if(ifr==1)
        legend(cl,'Location','northwest')
    end
    %
    subplot(2,2,3)
    xlabel('days')
    ylabel('numbers')
    title('linear-plot of numbers')
    axis([0 timex(end) 0 N/nzoom/10])
    if(ifr==1)
        legend(cl,'Location','northwest')
    end
    %
    subplot(2,2,4)
    xlabel('days')
    ylabel('numbers')
    title('linear-plot of numbers')
    axis([0 timex(end) 0 N/nzoom])
    if(ifr==1)
        legend(cl,'Location','northwest')
    end
    
    % surface plot of the local infection rate
    figure(12); surf(ic/Ncell)
end

    figure(11)
    cparas0 = [ ' mortality=' num2str(mort) ' days-sick=' num2str(sdays) ' action=' num2str(action1) ];
    cparas1 = [ ' cellsize=' num2str(idist) ' step=' num2str(step/4) ' p-fly=' num2str(probfly/dt) ' step-fly=' num2str(stepfly/4) ];
    h1z=1.15; % set first vertical position here
    h2z=1.09; % set second vertical position here
    h1x=-0.2; % set left-shift 
    h2x=1.2; % set left-shift 
    subplot(2,2,1)
    f1h0=text('string',cparas0,'interpreter','latex','fontsize',11,'units','norm','pos',[h1x h1z]);
    f1h1=text('string',cparas1,'interpreter','latex','fontsize',11,'units','norm','pos',[h2x h1z]);
    %title([cparas0 cparas1])
    cparas0
    cparas1
    ifrate/dt
    
    if(auto_mode==0)
        cmanual = input('Do you want to print? Y/N [N]:','s');
        if isempty(cmanual)
            cmanual='N';
        end 
    else
        cmanual='Y'
    end
    if( cmanual=='Y' || cmanual=='y' )
        for ifig=10:11
            figure(ifig); 
            orient(ifig,'landscape'); 
            print(ifig,['cor1_f' num2str(ifig) '_run' num2str(irun) '.pdf'],'-dpdf', '-fillpage')
        end
    end
end



