%   Script Paris.m
%
%   Bike hiring in Paris
%
% J. Lunze
% 2.5.2021
% for NCSAppl: 3.6.2021
%
echo off
clear
close all

%%  Graph with geographic positions of the hiring points
%
figure(1)
% 30 vertices in the following positions
N1=30;
xcoord=[4.5 2.0 4.0 8.6 1.4 12.3 13.5 11.5 9.6 1.8 12.9 5.5 12.3 0.9 10.3 5.5 7.3 10.8 10.2 7.8 2.2 1.8 3.5 4.2 4.2 9.5 11.8 7.3 10.5 12.5];
ycoord=[7.0 4.0 1.0 6.7 7.2  3.3  3.3  2.8 1.3 1.8  6.1 7.2  2.1 4.1  6.8 3.4 6.8  1.7  1.9 0.8 3.0 5.5 4.9 4.5 3.5 5.8  6.9 1.3  4.5  1.5];
%  Generation of the adjacency matrix of a random graph
p=0.04;
%Paris1=rng;                         % save s to restore the generator setting
%save('Paris1', 'Paris1');
load('Paris1');       % Seed (initialisation) of the number generator
rng(Paris1);
Adj=ceil(rand(N1,N1)-(1-p)*ones(N1,N1));
Adj=Adj-diag(diag(Adj));   % set the diagonal elements zero
G=digraph(Adj);
plot(G, 'XData', xcoord, 'YData', ycoord, 'LineWidth', 3, 'MarkerSize', 8,...
    'ArrowSize', 14, 'NodeFontName', 'Euclid', 'NodeFontSize', 14,...
    'NodeColor', 'black');
xleer;
yleer;
latextitle(['$${\tt Paris:}$$ Graph for $$p=$$', num2str(p)]);

%%  Second realisation
%
figure(2)
%Paris2=rng;                         % save s to restore the generator setting
%save('Paris2', 'Paris2');
load('Paris2');       % Seed (initialisation) of the number generator
rng(Paris2);
Adj=ceil(rand(N1,N1)-(1-p)*ones(N1,N1))';
       %transposition due to the MATLAB notation of adjacency matrices
Adj=Adj-diag(diag(Adj));   % set the diagonal elements zero
G=digraph(Adj);
plot(G, 'XData', xcoord, 'YData', ycoord, 'LineWidth', 3, 'MarkerSize', 8,...
    'ArrowSize', 14, 'NodeFontName', 'Euclid', 'NodeFontSize', 14,...
    'NodeColor', 'black');
xleer;
yleer;
latextitle(['$${\tt Paris:}$$ Graph for $$p=$$', num2str(p)]);

%%  Simulation of the bike hiring over one day
%
figure(3)
kend=34;
%Paris3=rng;                         % save s to restore the generator setting
%save('Paris3', 'Paris3');
load('Paris3');       % Seed (initialisation) of the number generator
rng(Paris3);
x0=6*ones(N1,1);
X=[];
X(1,:)=x0';
xold=X(1,:)';
for k1=1:kend
    Adj=ceil(rand(N1,N1)-(1-p)*ones(N1,N1));
    Adj=Adj-diag(diag(Adj));   % set the diagonal elements zero
    xnew=xold+Adj*ones(N1,1) - Adj'*ones(N1,1);
    X(k1+1,:)=xnew';
end
twosigma=2*sqrt((N1-1)*p*(1-p));
xmin=x0-twosigma;
xmax=x0+twosigma;
fill([0 kend+2 kend+2 0 0], [xmin xmin xmax xmax xmin], [0.9 0.9 0.9]);
hold on
plot([1:kend+1],X, 'b');
plot([0, kend+2], [x0, x0], 'k:');
latexxlabel('$$k$$');
latexylabel('$$x_i\;$$');
rotateY;
axis([0 kend+2 0 12]);
xticks([0 5 10 15 20 25 30 35]);
yticks([0 5 10 15 20]);
hold off















%%   Figures
epsfigc15('Paris');