function fig10_22

% y(t+dt)=Ax(t)

% measure dt in days
dt=7;

% https://usafacts.org/visualizations/coronavirus-covid-19-spread-map/

NY=importdata('covid_data.txt');
NY=NY';
[N M]=size(NY)

% XD = new cases; YD = new cases at t+dt
XD=NY(2:N-dt,:)-NY(1:N-dt-1,:);
YD=NY(dt+2:N,:)-NY(dt+1:N-1,:);
[n m]=size(XD)

%  data for table 10.5
% NY(1:4,1:3)
% XD(1:3,1:3)
% YD(n,1:3)
% XD(1:3,m-1:m)
% YD(n,m-1:m)

%%%% use normalized data
meansX=mean(XD);
meansY=mean(YD);
X=XD-meansX;
Y=YD-meansY;
for j=1:m
    XM(j)=norm(X(:,j),inf);
    YM(j)=norm(Y(:,j),inf);
    X(:,j)=X(:,j)/XM(j);
    Y(:,j)=Y(:,j)/YM(j);
end

% use the pseudo-inverse
% G=X'*X;
% cond_G=cond(G,inf)
% GI=inv(G);
% P=Y'*X*GI;

% use the SVD
[Uq,Sq,Vq]=svd(X);
SS=zeros(n,m);
for i=1:m
    SS(i,i)=1/Sq(i,i);
end
P=Y'*Uq*SS*Vq';

SSx=diag(1./XM);
Sy=diag(YM);
A=Sy*P*SSx;
b=meansY'-A*meansX';
%norm(b,inf)

%%% check on accuracy of y=Ax+b assumption
%%% Albany = 1
%%% Erie = 15
%%% Rensselaer = 42
%%% ulster = 56
%%% saratoga=46, king=24, westchester=60

r1=1; r2=15; r3=56;

figure(1)
clf
% get(gcf,'Position')
set(gcf,'Position', [4   925   560   420])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

subaxis(3,1,1,1,'MT',0.01,'MB',0.06,'MR',-0.001,'ML',0.06,'P',0.04)
r=r1;
yp=zeros(n,1);
for i=1:n
    yp(i)=dot(A(r,:),XD(i,:))+b(r);
end
max_yp=max(yp)
t=dt+1:n+dt;
plot(t,yp,'r','LineWidth',1.5)
hold on
plot(t,YD(:,r),'--b','LineWidth',1.5)
%xlabel('x-axis')
ylabel('Albany')
xlim([1 275])
ylim([0 300])
set(gca,'xTick',[1 91 181 274])
set(gca,'xTickLabel',{'April 1';'June 1';'Oct 1';'Jan 1'})
set(gca,'yTick',[0 150 300])
grid on
box on
set(gca,'FontSize',14,'FontWeight','bold')

%error1=norm(yp-YD(:,r)+XD(:,r),inf)

subaxis(3,1,1,2)
r=r2;
yp=zeros(n,1);
for i=1:n
    yp(i)=dot(A(r,:),XD(i,:))+b(r);
end
max_yp=max(yp)
plot(t,yp,'r','LineWidth',1.5)
hold on
plot(t,YD(:,r),'--b','LineWidth',1.5)
%xlabel('x-axis')
ylabel('Erie')
xlim([1 275])
ylim([0 800])
set(gca,'xTick',[1 91 181 274])
set(gca,'xTickLabel',{'April 1';'June 1';'Oct 1';'Jan 1'})
set(gca,'yTick',[0 400 800])
grid on
box on
legend({' Predicted',' Actual'},'Location','North')
set(gca,'FontSize',14,'FontWeight','bold')

%error2=norm(yp-YD(:,r)+XD(:,r),inf)


subaxis(3,1,1,3)
r=r3;
yp=zeros(n,1);
for i=1:n
    yp(i)=dot(A(r,:),XD(i,:))+b(r);
end
plot(t,yp,'r','LineWidth',1.5)
hold on
plot(t,YD(:,r),'--b','LineWidth',1.5)
xlabel('t-axis')
ylabel('Ulster')
xlim([1 275])
ylim([0 120])
set(gca,'yTick',[0 60 120])
set(gca,'xTick',[1 91 181 274])
set(gca,'xTickLabel',{'April 1';'June 1';'Oct 1';'Jan 1'})
grid on
box on
set(gca,'FontSize',14,'FontWeight','bold')

%error3=norm(yp-YD(:,r)+XD(:,r),inf)

%exportgraphics(gcf,'/Users/mark/Desktop/covid.eps')


figure(2)
clf
sings=svd(P);
fprintf('sigma(i)/sigma(1) = %6.2f, ',sings(1:4)/sings(1))
fprintf('\n\n')

% get(gcf)
set(gcf,'Position', [4 647 592 204])
plot(1:length(sings),sings/sings(1),'--o','LineWidth',1.5,'MarkerSize',8)
%axis([0 length(error) 0 0.1])
xlabel('k-axis')
ylabel('sigma/sigma(1)')
grid on
box on
set(gca,'FontSize',14,'FontWeight','bold')












