% Plot coronavirus cases for important areas

% remove old data
system('rm *.csv');

% the original data
system('cd /Users/fletch/Shared/src/covid-19-uk-data; git pull');
system('cp /Users/fletch/Shared/src/covid-19-uk-data/data/*.csv .');

% Lothian
system('grep Lothian covid-19-cases-uk.csv > lothian.csv');
lothian1 = readtable("lothian.csv");
lothian = table2array(lothian1(:,5:5))'; % extract cases
tsl = timeseries(lothian,1:length(lothian)); % TODO maybe get timeseries from input

tsl.Name = 'Lothian Cases';
tsl.TimeInfo.Units = 'days';
tsl.TimeInfo.StartDate = '01-March-2020';
tsl.TimeInfo.Format = 'dd mmm';

figure('Name', 'Lothian Coronavirus Cases')
figure(1)
plot(tsl)
title("Lothian Coronavirus Cases " + date)
saveas(gcf, 'lothian.png')

% Greater Glasgow & Clyde
system('grep Glasgow covid-19-cases-uk.csv > glasgow.csv');
ggc1 = readtable("glasgow.csv");
% ggc = str2double(table2array(ggc1(:,5:5))'); % extract cases allowing for
% thousands comma
ggc = table2array(ggc1(:,5:5))';
tsg = timeseries(ggc,1:length(ggc));

tsg.Name = 'Glasgow Cases';
tsg.TimeInfo.Units = 'days';
tsg.TimeInfo.StartDate = '01-March-2020';
tsg.TimeInfo.Format = 'dd mmm';

figure('Name', 'GG&C Coronavirus Cases')
figure(2)
plot(tsg)
title("GG&C Coronavirus Cases " + date)
saveas(gcf, 'glasgow.png')


% East Sussex
system('grep "East Sussex" covid-19-cases-uk.csv > esussex.csv');
esussex1 = readtable("esussex.csv");
esussex = table2array(esussex1(:,5:5))'; % extract cases
tss = timeseries(esussex,1:length(esussex));

tss.Name = 'East Sussex Cases';
tss.TimeInfo.Units = 'days';
tss.TimeInfo.StartDate = '05-March-2020';
tss.TimeInfo.Format = 'dd mmm';

figure('Name', 'East Sussex Coronavirus Cases')
figure(3)
plot(tss)
title("East Sussex Coronavirus Cases " + date)
saveas(gcf, 'east-sussex.png')

%
% Make graphs showing how quickly cases are doubling
%

lothDoub = covidDoubling(lothian);
figure('Name', 'Lothian Doubling')
lothDoub.Name = 'Lothian Days to Double Cases';
plot(lothDoub)
title("Lothian Days to Double Cases " + date)
saveas(gcf, 'lothian-days.png')

glasgowDoub = covidDoubling(ggc);
figure('Name', 'Glasgow Doubling')
glasgowDoub.Name = 'Glasgow Days to Double Cases';
plot(glasgowDoub)
title("GG&C Days to Double Cases " + date)
saveas(gcf, 'glasgow-days.png')


esussexDoub = covidDoubling(esussex);
figure('Name', 'East Sussex Doubling')
esussexDoub.Name = 'East Sussex Days to Double Cases';
esussexDoub.TimeInfo.StartDate = '05-March-2020';
plot(esussexDoub)
title("East Sussex Days to Double Cases " + date)
saveas(gcf, 'east-sussex-days.png')

% compare deaths between countries in UK
system('tail -1 covid-19-totals-england.csv > englast.csv');
englast = csvread("englast.csv",0,1);
engpop = 55.93; % in millions from http://www.ukpopulation.org/england-population/
engdpm = englast(1,3)/engpop; % deaths per million population

system('tail -1 covid-19-totals-northern-ireland.csv > nilast.csv');
nilast = csvread("nilast.csv",0,1);
nipop = 1.908250; % from http://www.ukpopulation.org/
nidpm = nilast(1,3)/nipop; % deaths per million population

system('tail -1 covid-19-totals-scotland.csv > scotlast.csv');
scotlast = csvread("scotlast.csv",0,1);
scotpop = 5.50; % from http://www.ukpopulation.org/
scotdpm = scotlast(1,3)/scotpop; % deaths per million population

system('tail -1 covid-19-totals-wales.csv > wallast.csv');
wallast = csvread("wallast.csv",0,1);
walpop = 3.230490; % from http://www.ukpopulation.org/
waldpm = wallast(1,3)/walpop; % deaths per million population

figure('Name', 'UK Countries Deaths Per Million Population')
title('UK Countries Deaths Per Million Population')
countries = categorical({'England', 'NI', 'Scotland', 'Wales'});
bar(countries, [engdpm nidpm scotdpm waldpm])
title("UK Countries Deaths Per Million Population " + date)
saveas(gcf, 'ukdeaths.png')

msgbox('remember to copy PNGs to the web site', 'Cases', 'warn')

%
% FUNCTIONS
% 

function dts = covidDoubling(casests)
    % Transform cases TimeSeries to show doubling rate
    %
    % iterate thro the TimeSeries casests
    %
    for i = 1:length(casests)
    % if index is < 5 then 0
    if i < 5
        doubling = 0;
        % if value < 100 then 0  // ignore early values
    elseif casests(i) < 100
        doubling = 0;
        %
        % if the five previous values are less than the current value 
        %   then the values are increasing
        % else if the five previous values are more than the current value
        %   then the values are decreasing 
        %   (and values will be shown as negative)
        % else treat as increasing
    else
        incdec = 1; % incdec is 1 if increasing, -1 if decreasing
        for j = i-6:i-1
            if casests(i) < casests(j)
                incdec = -1;
            end
        end

        %
        % look backwards through the TimeSeries and count the number of 
        % days since the value was half (or double) the current value.
        if incdec == 1
            halfCurrent = casests(i) / 2; % increasing
        else
            halfCurrent = casests(i) * 2; % decreasing
        end
        doubling = 0;
        k = i;
        while casests(k) > halfCurrent
            doubling = doubling + 1;
            k = k - 1;
        end

        % Store the count as positive if increasing and negative if decreasing
        doubling = doubling * incdec;
    end
    doublings(i) = doubling;
    end
    dts = timeseries(doublings, 1:length(doublings));
    % dts.Name = 'Lothian Days to Double Cases';
    dts.TimeInfo.Units = 'days';
    dts.TimeInfo.StartDate = '01-March-2020';
    dts.TimeInfo.Format = 'dd mmm';
end