path="c:/users/dre2/desktop/DRYGINresubmission/";

(* Build a new data set by resampling the data under the assumption of zero epistasis and t distributed standard deviations. Group D *)

singlekor=Import[StringJoin[path,"FILES2012/a_tsq.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_genuinearray.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoD.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaD.dat"]];

fakedoubleko=Table["NaN",{Length[singlekor]},{Length[singlekoc]}];

For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],j++,
If[doubleko[[i,j]]=="NaN"||doublekosigma[[i,j]]=="NaN",Continue[]];
If[doubleko[[i,j]]<0,Continue[]];
sigma=doublekosigma[[i,j]];
If[sigma<(10^(-6)),sigma=(10^(-6))];
fakedoubleko[[i,j]]=RandomReal[StudentTDistribution[singlekor[[i,2]] singlekoc[[j,2]],(sigma/2),3]];
];
]
Export[StringJoin[path,"FILES2012/j4b_doublekoD_fake_studentt.dat"],fakedoubleko]

(* Create the epistasis matrix for group D: tsq x genuinearray *)

singlekor=Import[StringJoin[path,"FILES2012/h4a_tsq_fake_studentt.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/h4a_genuinearray_fake_studentt.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/j4b_doublekoD_fake_studentt.dat"]];

eps=Table[0,{Length[singlekor]},{Length[singlekoc]}];

For[n1=1,n1≤Length[singlekor],n1++,
For[n2=1,n2≤Length[singlekoc],n2++,

a=doubleko[[n1,n2]];

If[a=="NaN",
eps[[n1,n2]]="NaN",
If[a≤0,
Print["Non positive growthrate for the pair (",n1,",",n2,"): ",doubleko[[n1,n2]]],
eps[[n1,n2]]=Log[2,(a/(singlekor[[n1,2]]singlekoc[[n2,2]]))];
];
];

];
]
Export[StringJoin[path,"FILES2012/j4c_epsD_fake_studentt.dat"],eps]

(* expbinning[r] implements the exponential binning. The value returned, is an integer from 1 to dim=17. *)

expbinning[r_]:=Module[{bin},

bin=0;

If[r<-(2^1),bin=1];
If[r≥ -(2^1)&&r<-(2^(-0)),bin=2;Return[bin];];
If[r≥ -(2^(-0))&&r<-(2^(-1)),bin=3;Return[bin];];
If[r≥ -(2^(-1))&&r<-(2^(-2)),bin=4;Return[bin];];
If[r≥ -(2^(-2))&&r<-(2^(-3)),bin=5;Return[bin];];
If[r≥ -(2^(-3))&&r<-(2^(-4)),bin=6;Return[bin];];
If[r≥ -(2^(-4))&&r<-(2^(-5)),bin=7;Return[bin];];
If[r≥ -(2^(-5))&&r<-(2^(-6)),bin=8;Return[bin];];
If[r≥ -(2^(-6))&&r<-(2^(-7)),bin=9;Return[bin];];
If[r≥ -(2^(-7))&&r<-(2^(-8)),bin=10;Return[bin];];
If[r≥ -(2^(-8))&&r<(2^(-8)),bin=11;Return[bin];];
If[r≥ (2^(-8))&&r<(2^(-7)),bin=12;Return[bin];];
If[r≥ (2^(-7))&&r<(2^(-6)),bin=13;Return[bin];];
If[r≥ (2^(-6))&&r<(2^(-5)),bin=14;Return[bin];];
If[r≥ (2^(-5))&&r<(2^(-4)),bin=15;Return[bin];];
If[r≥ (2^(-4))&&r<(2^(-3)),bin=16;Return[bin];];
If[r≥ (2^(-3)),bin=17;Return[bin];];

bin
]

dim=17;

(* Bin each epistatic interaction eps(f1,f2) as a function of df1 and df2. For each bin compute the number of points in each bin and the first moment and the square of the second moment. Group D. *)

nnanpoints=Table[0,{dim},{dim}];
npoints=Table[0,{dim},{dim}];
delta1=Table[0,{dim},{dim}];
delta2=Table[0,{dim},{dim}];
binx=Table[0,{dim}];
nx=Table[0,{dim}];

q25=Table[0,{dim},{dim}];
q50=Table[0,{dim},{dim}];
q75=Table[0,{dim},{dim}];
q90=Table[0,{dim},{dim}];

For[i=1,i≤dim,i++,
For[j=1,j≤dim,j++,

quantilelist={};

(* Group D = tsq x genuinearray. *)

singlekor=Import[StringJoin[path,"FILES2012/h4a_tsq_fake_studentt.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/h4a_genuinearray_fake_studentt.dat"]];
eps=Import[StringJoin[path,"FILES2012/j4c_epsD_fake_studentt.dat"]];

(* Create look-up tables associating to each mutant its log-growth rate and its bin. *)

logsinglekor=Table[0,{Length[singlekor]}];
binsinglekor=Table[0,{Length[singlekor]}];
For[n=1,n≤Length[singlekor],n++,
logsinglekor[[n]]=Log[2,singlekor[[n,2]]];
binsinglekor[[n]]=expbinning[logsinglekor[[n]]];
];

logsinglekoc=Table[0,{Length[singlekoc]}];
binsinglekoc=Table[0,{Length[singlekoc]}];
For[n=1,n≤Length[singlekoc],n++,
logsinglekoc[[n]]=Log[2,singlekoc[[n,2]]];
binsinglekoc[[n]]=expbinning[logsinglekoc[[n]]];
];

For[n1=1,n1≤Length[singlekor],n1++,
For[n2=1,n2≤Length[singlekoc],n2++,

df1=logsinglekor[[n1]];
df2=logsinglekoc[[n2]];
a=binsinglekor[[n1]];
b=binsinglekoc[[n2]];

(*If[a≠i||b≠j,Continue[]];*)
If[!((a==i&&b==j)||(a==j&&b==i)),Continue[]];

If[eps[[n1,n2]]=="NaN",nnanpoints[[a,b]]++;nnanpoints[[b,a]]++;Continue[]];

(* each double-knockout is counted as two points *)
npoints[[a,b]]++;
delta1[[a,b]]+=eps[[n1,n2]];
delta2[[a,b]]+=((eps[[n1,n2]])^2);

npoints[[b,a]]++;   (* different from A to double count correctly *)
delta1[[b,a]]+=eps[[n1,n2]];
delta2[[b,a]]+=((eps[[n1,n2]])^2);

AppendTo[quantilelist,((eps[[n1,n2]])^2)];

binx[[a]]+= df1;
nx[[a]]++;
binx[[b]]+= df2;
nx[[b]]++;
];
];

If[Length[quantilelist]>1,
q25[[i,j]]=Quantile[quantilelist,1/4];
q50[[i,j]]=Quantile[quantilelist,1/2];
q75[[i,j]]=Quantile[quantilelist,3/4];
q90[[i,j]]=Quantile[quantilelist,9/10];
q25[[j,i]]=q25[[i,j]];
q50[[j,i]]=q50[[i,j]];
q75[[j,i]]=q75[[i,j]];
q90[[j,i]]=q90[[i,j]];
];

Print[i,"  ",j];
];
];

Export[StringJoin[path,"FILES2012/j4d_nnanpointsD_fake_studentt.dat"],nnanpoints];
Export[StringJoin[path,"FILES2012/j4d_npointsD_fake_studentt.dat"],npoints];
Export[StringJoin[path,"FILES2012/j4d_delta1D_fake_studentt.dat"],delta1];
Export[StringJoin[path,"FILES2012/j4d_delta2D_fake_studentt.dat"],delta2];
Export[StringJoin[path,"FILES2012/j4d_binxD_fake_studentt.dat"],binx];
Export[StringJoin[path,"FILES2012/j4d_nxD_fake_studentt.dat"],nx];

Export[StringJoin[path,"FILES2012/j4d_q25D_fake_studentt.dat"],q25];
Export[StringJoin[path,"FILES2012/j4d_q50D_fake_studentt.dat"],q50];
Export[StringJoin[path,"FILES2012/j4d_q75D_fake_studentt.dat"],q75];
Export[StringJoin[path,"FILES2012/j4d_q90D_fake_studentt.dat"],q90];

(* Go back from exponential binning to linear scale and plot: number of points, average deviation from the linear model, variance and standard deviation. Group D *)

nnanpoints=Import[StringJoin[path,"FILES2012/j4d_nnanpointsD_fake_studentt.dat"]];
npoints=Import[StringJoin[path,"FILES2012/j4d_npointsD_fake_studentt.dat"]];
delta1=Import[StringJoin[path,"FILES2012/j4d_delta1D_fake_studentt.dat"]];
delta2=Import[StringJoin[path,"FILES2012/j4d_delta2D_fake_studentt.dat"]];
binx=Flatten[Import[StringJoin[path,"FILES2012/j4d_binxD_fake_studentt.dat"]]];
nx=Flatten[Import[StringJoin[path,"FILES2012/j4d_nxD_fake_studentt.dat"]]];

q25=Import[StringJoin[path,"FILES2012/j4d_q25D_fake_studentt.dat"]];
q50=Import[StringJoin[path,"FILES2012/j4d_q50D_fake_studentt.dat"]];
q75=Import[StringJoin[path,"FILES2012/j4d_q75D_fake_studentt.dat"]];
q90=Import[StringJoin[path,"FILES2012/j4d_q90D_fake_studentt.dat"]];

meanx=Table[0,{dim}];
For[i=1,i≤dim,i++,
If[nx[[i]]≠0,meanx[[i]]=((binx[[i]])/(nx[[i]])),meanx[[i]]=Null];
]

nangrid={};
eps0={};
eps1={};
eps2={};
epssigma={};

eps25={};
eps50={};
eps75={};
eps90={};

For[i=1,i≤dim,i++,
For[j=1,j≤dim,j++,

x=meanx[[i]];
y=meanx[[j]];

If[npoints[[i,j]]>4,   (* Notice the 4 here *)
AppendTo[nangrid,{x,y,Log[10,nnanpoints[[i,j]]]}];
AppendTo[eps0,{x,y,Log[10,npoints[[i,j]]]}];
AppendTo[eps1,{x,y,((delta1[[i,j]])/(npoints[[i,j]]))}];
AppendTo[eps2,{x,y,((delta2[[i,j]])/(npoints[[i,j]]))}];
AppendTo[epssigma,{x,y,Sqrt[((delta2[[i,j]])/(npoints[[i,j]]))]}];

AppendTo[eps25,{x,y,Sqrt[q25[[i,j]]]}];
AppendTo[eps50,{x,y,Sqrt[q50[[i,j]]]}];
AppendTo[eps75,{x,y,Sqrt[q75[[i,j]]]}];
AppendTo[eps90,{x,y,Sqrt[q90[[i,j]]]}];
];
];
];

(*ListPlot3D[nangrid,PlotRange->All,AxesStyle->18]*)
ListPlot3D[eps0,PlotRange->All,AxesStyle->18]
ListPlot3D[eps1,PlotRange->All,AxesStyle->18]
ListPlot3D[eps2,PlotRange->All,AxesStyle->18]
prms=ListPlot3D[epssigma,PlotRange->All,AxesStyle->18,PlotStyle->Green]

p25=ListPlot3D[eps25,PlotRange->All,AxesStyle->18,PlotStyle->Blue]
p50=ListPlot3D[eps50,PlotRange->All,AxesStyle->18,PlotStyle->Purple]
p75=ListPlot3D[eps75,PlotRange->All,AxesStyle->18,PlotStyle->Blue]
p90=ListPlot3D[eps90,PlotRange->All,AxesStyle->18,PlotStyle->Cyan]

Show[p25,p50,p75,p90,prms]