%Code developed by 
%Amanda M. Wilson; Marco-Felipe King, PhD; and Martin Lopez-Garcia, PhD
%through University of Arizona and University of Leeds collaboration

%This code includes

%1 - Entering experimental data
%2 - Approximate Bayesian computation approach (ABC)
%3 - QMRA application
%4 - Code for Figure 1

%%
clear all
%% 1 - Entering experimental data

%col 1: concentration on the hand before the contact
%col 2: concentration on the surface before the contact
%col 3: concentration on the hand after the contact
%col 4: concentration on the surface after the contact

%These data are also available in a separate file,
%titled, "experimental_data.xlsx"

C=[
8.25E+05,	2.40E+05,	6.67E+05,   2.12E+05 
1.25E+05,	5.35E+06,	3.32E+05,   4.30E+06
7.30E+05,	2.72E+04,	6.51E+05,   2.71E+04 
5.20E+04,	3.10E+06,	4.93E+04,   3.05E+06 
2.34E+05,	3.50E+06,	2.39E+05,   9.16E+05
8.25E+04,	2.50E+06,	5.92E+04,   2.90E+06
7.70E+05,	5.30E+07,	9.28E+05,   5.35E+07 
3.05E+05,	1.41E+08,	1.67E+06,   1.71E+08 
1.83E+02,	3.80E+06,	4.92E+04,   2.15E+06 
2.70E+02,	2.30E+06,	1.12E+04,   1.46E+06 
3.70E+07,	3.35E+06,	1.59E+07,   1.98E+06 
5.20E+05,	7.35E+02,	5.75E+05,   4.52E+03 
2.52E+07,	1.75E+06,	2.67E+07,   2.14E+06 
2.02E+08,	1.32E+06,	2.07E+08,   2.54E+06 
3.35E+05,	4.50E+02,	4.45E+05,   3.81E+03 
4.80E+05,	5.00E+00,	6.40E+05,   1.15E+03 
];

%% 2 - Approximate Bayesian computation approach (ABC)

%Number of iterations
n=1e07;

%---------------------Prior Distributions----------------------------------

%Prior distribution for lambda
pd=makedist('uniform',0.0001,0.406);
lambdaU=random(pd,n,1);

%Prior distributions for fingertip swabbing efficiency
pd=makedist('uniform',0,1);
sEffhand=random(pd,n,1);
%and for surface swabbing efficiency
sEffsurf=random(pd,n,1);

%---Estimating fingertip and surface viral concentrations after contact---

%Calculating estimated concentrations on finger tip and surface using
Calc_C_H=((C(:,1)./sEffhand'-(lambdaU'.*(C(:,1)./sEffhand'-C(:,2)./sEffsurf'))));
Calc_C_S=((C(:,2)./sEffsurf'-(lambdaU'.*(C(:,2)./sEffsurf'-C(:,1)./sEffhand'))));

%-------------------Calculating Euclidean distance-------------------------

%calculating individual Euclidean distances for fingertip and surface
tmpH=log10(Calc_C_H(1:end,:).*sEffhand')-log10(C(:,3));
tmpS=log10(Calc_C_S(1:end,:).*sEffsurf')-log10(C(:,4));

%Combined Euclidean distance that takes fingerip and surface into account
dist= sqrt(sum(tmpH.^2+tmpS.^2)); %Combined

%-------------------Posterior distributions--------------------------------

%sorting distances and subsetting to get posterior distributions

%restriction
restrict=0.001;

%indices of sorted data
[~, Sort]=sort(dist);

%indices of best distances, based on our restriction level
[best]=Sort(1:(restrict*n));

%--------subsetting data - combinatins with best distances-----------------

%estimated concentration on fingertip after contact
Calc_C_HSorted=Calc_C_H(:,best); 

%estimated concentration on surface after contact
Calc_C_SSorted=Calc_C_S(:,best);

%posterior distribution for lambda
lambdaUSorted=lambdaU(best);

%posterior distribution for fingertip swabbing efficiency
sEffhandSorted=sEffhand(best);

%posterior distribution for surface swabbing efficiency
sEffsurfSorted=sEffsurf(best);

%estimated concentration on hand accounting for swabbing efficiency
Calc_C_HSorted_swab=Calc_C_HSorted.*sEffhandSorted';

%estimated concentration on surface accounting for swabbing efficiency
Calc_C_SSorted_swab=Calc_C_SSorted.*sEffsurfSorted';

%writing csv files for ABC step
csvwrite('Calc_C_HSorted.csv',Calc_C_HSorted_swab);
csvwrite('Calc_C_SSorted.csv',Calc_C_SSorted_swab);
csvwrite('sEffhandSorted.csv',sEffhandSorted);
csvwrite('sEffsurfSorted.csv',sEffsurfSorted);
csvwrite('lambda.csv',lambdaUSorted);
csvwrite('C_data.csv',C);

%% 3 - QMRA application

%-----------------Input parameter Distributions----------------------------

%hand-to-mouth transfer efficiency
pd=makedist('normal',.3390,0.1);
mouthte=random(truncate(pd,0,1),n*restrict,1);

%concentration on the fingertip before the contact
pd=makedist('uniform',-2,2);
H0=10.^random(pd,n*restrict,1);

%concentration on the surface before the ocntact
pd=makedist('uniform',-2,2);
S0=10.^random(pd,n*restrict,1);

%----Singe-fomite and -mouth contact exposure & infection risk estimate----

%estimating conc on hand after surface contact using posterior lambda
%distribution
Hpredict=H0-(lambdaUSorted.*(H0-S0));

%predicting dose assuming whole finger tip used in surface contact is also
%made in hand-to-mouth contact

%Predicted doses here will be related to infection risk in .r code
Dosepredict=mouthte.*Hpredict;
%----------------------saving output---------------------------------------

%writing csv files for QMRA step
csvwrite('Dosepredict.csv',Dosepredict);
csvwrite('Hpredict.csv',Hpredict);
csvwrite('S0.csv',S0);
csvwrite('H0.csv',H0);
csvwrite('mouthte.csv',mouthte);
csvwrite('lambda.csv',lambdaUSorted);
%% 4 - Code for Figure 1

varplot=3;
matrix_for_plot=zeros(n*restrict,varplot);
matrix_for_plot(:,1)=sEffhandSorted;
matrix_for_plot(:,2)=sEffsurfSorted;
matrix_for_plot(:,3)=lambdaUSorted;

figure;
plotmatrix(matrix_for_plot);
[~,ax]=plotmatrix(matrix_for_plot);
ax(1,1).YLabel.String='Fingertip Swabbing Efficiency'; 
ax(2,1).YLabel.String='Surface Swabbing Efficiency'; 
ax(3,1).YLabel.String='{\lambda}'; 
ax(end,1).XLabel.String='Fingertip Swabbing Efficiency'; 
ax(end,2).XLabel.String='Surface Swabbing Efficiency'; 
ax(end,3).XLabel.String='{\lambda}';  
