function fig10_23

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

% measure dt in days
dt=7;

% COVID data
% https://usafacts.org/visualizations/coronavirus-covid-19-spread-map/
% open Download Data/Known Cases and save as a text file covid_data
T = readtable('covid_data3.txt');
[Nd,Md]=size(T);

% NY rows 1865 to 1926
% April 1, 2020 = column 75

% T(1865:1926,75:80)
%
% pause

NY=table2array(T(1865:1926,75:Md))';
[N M]=size(NY)

n=N-dt;
m=M;
XD=NY(1:n,:);
YD=NY(dt+1:n+dt,:);

%return

idata=1;
if idata==1
    %%%% use normalized data
    X=zeros(n,m);
    Y=zeros(n,m);
    meansX=mean(XD);
    meansY=mean(YD);
    for j=1:m
        X(:,j)=XD(:,j)-meansX(j)*ones(n,1);
        Y(:,j)=YD(:,j)-meansY(j)*ones(n,1);
    end
    for j=1:m
        XM(j)=max(abs(X(:,j)));
        YM(j)=max(abs(Y(:,j)));
        X(:,j)=X(:,j)/XM(j);
        Y(:,j)=Y(:,j)/YM(j);
    end
elseif idata==2
    % use unnormalized data
    X=XD;
    Y=YD;
    meansX=zeros(1,m);
    meansY=meansX;
    XM=ones(1,n);
    YM=ones(1,n);
end

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

%%% Albany = 1
%%% Erie = 15
%%% Rensselaer = 42
%%% ulster = 56
%%% saratoga=46, king=24, westchester=60
r=1

P(r,r)
ss=sort(P(r,:))
for ia=1:length(ss)
    if P(r,ia) == ss(1)
        [ia P(r,ia)]
    elseif P(r,ia) == ss(length(ss))
        [ia P(r,ia)]
    end
end


% NY state GIS data
% http://gis.ny.gov/gisdata/inventories/details.cfm?DSID=927
% Download NYS Civil Boundaries (Shape)
%data = shaperead('Counties.shp');
data = shaperead('Counties_Shoreline');
%mapshow(data)

data2=data;
data2(23,:)=data(62,:);
data2(24:26,:)=data(23:25,:);
data2(27,:)=data(60,:);
data2(28:32,:)=data(26:30,:);
data2(33,:)=data(61,:);
data2(34:44,:)=data(31:41,:);
data2(45,:)=data(47,:);
data2(46:50,:)=data(42:46,:);
data2(51:62,:)=data(48:59,:);

data2(30,:)
data2(52,:)

clf
% get(gcf,'Position')
set(gcf,'Position', [3   925   560   420])
subaxis(1,1,1,1,'MT',-0.01,'MB',-0.02,'MR',-0.1,'ML',-0.1,'P',0.04)

pm=min(P(r,:))
pM=max(P(r,:))
for ic = 1:62
    %     ic
    %     data2(ic,:)
    if P(r,ic)>0
        gg=1-P(r,ic)/pM;  bb=gg;
        if ic==r
            mapshow(data2(ic,:),'FaceColor',[1 gg bb],'LineWidth',4);
        else
            mapshow(data2(ic,:),'FaceColor',[1 gg bb],'LineWidth',1);
        end
    else
        rr=1-P(r,ic)/pm;
        gg=rr;
        if ic==r
            mapshow(data2(ic,:),'FaceColor',[rr gg 1],'LineWidth',4);
        else
            mapshow(data2(ic,:),'FaceColor',[rr gg 1],'LineWidth',1);
        end
    end
    %pause
end

axis image off manual
xl = xlim
yl = ylim
xl(1)
say=['d = ', num2str(dt,'%i')];
text(xl(1)+0.1*(xl(2)-xl(1)),yl(1)+0.9*(yl(2)-yl(1)),say,'FontSize',30,'FontWeight','bold')












