Datadir= $HomeDirectory<>"\\Desktop\\AfterDisert.Work\\Chr3Potential\\";
SNPsFileNames=FileNames["*", Datadir];
EnvParDA={{FIN,{Alt,59.35},{TEM,3.604},{Ran,48.206},{Bac,21.435},{Vir,16.197},{Pro,7.775},{Hel,27.065},{Chi,2.112},{Pri,0.6},{Rod,11.509},{Sor,5.099},{UVB,"NA"},{WS2M,2.497},{RH2M,164.595},{PS,187.776}},{CHB,{Alt,134.56},{TEM,14.138},{Ran,74.709},{Bac,21.145},{Vir,15.612},{Pro,8.391},{Hel,24.658},{Chi,3.466},{Pri,0.884},{Rod,9.992},{Sor,1.360},{UVB,1958.224},{WS2M,2.969},{RH2M,66.192},{PS,99.556}},{LWK,{Alt,1391.19},{TEM,20.310},{Ran,122.904},{Bac,23.8},{Vir,21},{Pro,10.4},{Hel,26.20},{Chi,8.3},{Pri,4.5},{Rod,5.5},{Sor,0.6},{UVB,5977.528},{WS2M,1.607},{RH2M,68.97},{PS,85.22}},{TSI,{Alt,48.45},{TEM,12.728},{Ran,65.979},{Bac,23.529},{Vir,19.294},{Pro,9.412},{Hel,25.6},{Chi,8.24},{Pri,0.6},{Rod,13.34},{Sor,4.6},{UVB,539.326},{WS2M,2.017},{RH2M,74.610},{PS,97.421}},{JPT,{Alt,58.46},{TEM,13.035},{Ran,126.549},{Bac,22.538},{Vir,14.725},{Pro,8.963},{Hel,25.880},{Chi,2.648},{Pri,1.679},{Rod,8.151},{Sor,1.9},{UVB,1168.539},{WS2M,2.403},{RH2M,79.082},{PS,99.558}},{IBS,{Alt,441.462},{TEM,14.555},{Ran,43.058},{Bac,25.904},{Vir,20.239},{Pro,10.206},{Hel,28.246},{Chi,8.705},{Pri,.6},{Rod,7.924},{Sor,2.751},{UVB,1009.971},{WS2M,2.483},{RH2M,61.26},{PS,95.28}},{KHV,{Alt,5.373},{TEM,26.253},{Ran,158.593},{Bac,21.77},{Vir,17.354},{Pro,9.663},{Hel,16.188},{Chi,14.439},{Pri,3.101},{Rod,7.611},{Sor,1.648},{UVB,3953.286},{WS2M,2.292},{RH2M,76.658},{PS,100.177}},{MSL,{Alt,74.3},{TEM,25.805},{Ran,246.568},{Bac,21.741},{Vir,20.744},{Pro,7.440},{Hel,25.035},{Chi,10.1},{Pri,3.4},{Rod,9.8},{Sor,0.6},{UVB,4719.101},{WS2M,1.959},{RH2M,81.355},{PS,100.401}},{CLM,{Alt,1915.529},{TEM,17.569},{Ran,115.518},{Bac,18.99},{Vir,15.196},{Pro,7.322},{Hel,9.034},{Chi,21.118},{Pri,2.337},{Rod,6.015},{Sor,0.6},{UVB,5036.091},{WS2M,1.173},{RH2M,78.997},{PS,83.530}},{YRI,{Alt,79.60},{TEM,27.074},{Ran,122.873},{Bac,20.399},{Vir,17.673},{Pro,10.295},{Hel,19.599},{Chi,8.572},{Pri,3.662},{Rod,7.968},{Sor,0.6},{UVB,5115.505},{WS2M,1.670},{RH2M,83.621},{PS,99.997}},{PEL,{Alt,292.234},{TEM,18.393},{Ran,16.018},{Bac,6.215},{Vir,10.740},{Pro,6.410},{Hel,6.252},{Chi,8.632},{Pri,0.726},{Rod,2.296},{Sor,0.6},{UVB,7747.893},{WS2M,3.715},{RH2M,80.339},{PS,95.564}}};

Print["Defining data retrieval Functions from directory " ,Datadir];
AllelesSwitcher[alleles_,lst_]:=Which[
ToString[lst[[1,1]]]==ToString[alleles],lst,
True, ReplacePart[lst,{{1,1}->alleles,{2,1}->alleles,{1,2,3}->lst[[1,2,2]],{2,2,3}->lst[[2,2,2]],{1,2,2}->lst[[1,2,3]],{2,2,2}->lst[[2,2,3]]}]
];
RedOrBlue[dataset_,j_]:=Which[
((Length[dataset[[j,3,1,2]]]==3)||(dataset[[j,3,1,2]][[4]]==0.)),Blue,
True,Red
];
SNPdata[rsno_]:=Block[
{SNPfile=FileNames["*"<>ToString[rsno],Datadir][[1]]},
Get[SNPfile]
];
orderedAlleleData[dataset_]:=Block[{finalorderedalleledata={dataset[[1]]},
jfoad=2},
While[jfoad≤Length[dataset],
newpop={dataset[[jfoad,1]],dataset[[jfoad,2]]};
switchpop=AppendTo[newpop,AllelesSwitcher[finalorderedalleledata[[1,3,1,1]],dataset[[jfoad,3]]]];
AppendTo[finalorderedalleledata,switchpop];
jfoad++
];
finalorderedalleledata
];

AdjParametersAndPotentials[paramno_,SNPfileno_]:=Block[{populations=Transpose[EnvParDA][[1]],
paramsandnames=Transpose[EnvParDA][[1+paramno]],paramname=Transpose[EnvParDA][[1+paramno]][[1,1]],params=Transpose[Transpose[EnvParDA][[1+paramno]]][[2]],OAData=Which[
(dat=Get[SNPsFileNames[[SNPfileno]]])≠{},orderedAlleleData[dat],
True,{}
]},
(*This will be a loop over all populations*)
finaldata={};
For[popno=1,popno≤Length[populations],popno++,
If[((pos=Position[OAData,ToString[populations[[popno]]]])!={})&& (params[[popno]]≠"NA") ,
lstno=pos[[1,1]];
AppendTo[finaldata,{{populations[[popno]],RedOrBlue[OAData,lstno]},{params[[popno]],OAData[[lstno,3,1]]}}]
];
];
Sortedfinaldata=
Sort[finaldata,#1[[2,1]]≤#2[[2,1]] &];
SNPNum=Which[
OAData≠{},OAData[[1,1,3]],
True,StringDrop[SNPsFileNames[[SNPfileno]],StringPosition[SNPsFileNames[[SNPfileno]],"\\"][[-1,-1]]]
];
{{SNPNum,paramname},Sortedfinaldata}
];
SNPAndAllelicMusvsParam[paramno_,SNPfileno_]:=
Block[{paramsandpots=AdjParametersAndPotentials[paramno,SNPfileno],SNPPoints=Table[{AdjParametersAndPotentials[paramno,SNPfileno][[2,j,2,1]],AdjParametersAndPotentials[paramno,SNPfileno][[2,j,2,2,2,1]]},{j,1,Length[AdjParametersAndPotentials[paramno,SNPfileno][[2]]]}],Allele1Points=Table[{AdjParametersAndPotentials[paramno,SNPfileno][[2,j,2,1]],AdjParametersAndPotentials[paramno,SNPfileno][[2,j,2,2,2,2]]},{j,1,Length[AdjParametersAndPotentials[paramno,SNPfileno][[2]]]}],Allele2Points=Table[{AdjParametersAndPotentials[paramno,SNPfileno][[2,j,2,1]],AdjParametersAndPotentials[paramno,SNPfileno][[2,j,2,2,2,3]]},{j,1,Length[AdjParametersAndPotentials[paramno,SNPfileno][[2]]]}],Alleles=AdjParametersAndPotentials[paramno,SNPfileno][[2,1,2,2,1]],PopsAndColor=Transpose[AdjParametersAndPotentials[paramno,SNPfileno][[2]] ][[1]],SNPName=AdjParametersAndPotentials[paramno,SNPfileno][[1,1]],envpar=AdjParametersAndPotentials[paramno,SNPfileno][[1,2]]},
{{SNPName,envpar},PopsAndColor,SNPPoints,Alleles,Allele1Points,Allele2Points}
];
RMS[points_,Functions_]:=Sqrt[Sum[(points[[j,2]]-Functions[points[[j,1]]])^2/Length[points],{j,1,Length[points]}]];
muafunc[pa1_,Tp_,Mu_]:=Mu- Tp Log[2,pa1];
musfunc[pa1_,Tp_,Mu_]:=Mu- Tp(pa1 Log[2,pa1] + (1- pa1) Log[2,1-pa1] );
MyFun[type_][ptlist_,funcform_,funcparams_,var_]:=Block[{Potentialrange=Abs[Max[Transpose[ptlist][[2]]]-Min[Transpose[ptlist][[2]]]],ans=FindFit[ptlist,funcform,funcparams,var]},
sss[type]=funcform/.ans;
Ftfunc[type][xx_]:=sss[type]/.var->xx;
RMS[ptlist,Ftfunc[type]]/Potentialrange
];
SNPfuncPlot[Plotinfo_,func_,flagval_]:=Block[
{height=(Max[Transpose[Plotinfo[[3]]][[2]]]-Min[Transpose[Plotinfo[[3]]][[2]]]),
width=(Max[Transpose[Plotinfo[[3]]][[1]]]-Min[Transpose[Plotinfo[[3]]][[1]]])},
Plot[func[var],{var,Min[Transpose[Plotinfo[[3]]][[1]]],Max[Transpose[Plotinfo[[3]]][[1]]]},
PlotStyle->AbsoluteThickness[Min[{(Maxflagvalue/flagval),10}]],PlotRange->{{Min[Transpose[Plotinfo[[3]]][[1]]]-0.05width,Max[Transpose[Plotinfo[[3]]][[1]]]+0.05width},{(Min[Transpose[Plotinfo[[3]]][[2]]]-0.03*height),(Max[Transpose[Plotinfo[[3]]][[2]]]+0.05*height)}},AspectRatio->Full,AxesLabel->{Plotinfo[[1,2]],"μ"<>ToString[Plotinfo[[1,1]]]},PlotLabel->ToString[Plotinfo[[1,1]]]<>" Relative RMS ="<>ToString[N[Round[flagval*1000]/1000]],Epilog->Table[{Transpose[Plotinfo[[2]]][[2,i]],PointSize[0.01],Point[Plotinfo[[3,i]]],Text[Plotinfo[[2,i,1]],Plotinfo[[3,i]]+{0.0,0.03*height}]},{i,1,Length[Plotinfo[[2]]]}]]
];
MuAllele1FuncPlot[Plotinfo_,func_,flagval_]:=Block[
{height=(Max[Transpose[Plotinfo[[5]]][[2]]]-Min[Transpose[Plotinfo[[5]]][[2]]]),
width=(Max[Transpose[Plotinfo[[5]]][[1]]]-Min[Transpose[Plotinfo[[5]]][[1]]])},
Plot[func[var],{var,Min[Transpose[Plotinfo[[5]]][[1]]],Max[Transpose[Plotinfo[[5]]][[1]]]},
PlotStyle->AbsoluteThickness[Min[{(Maxflagvalue/flagval),10}]],PlotRange->{{Min[Transpose[Plotinfo[[5]]][[1]]]-0.05width,Max[Transpose[Plotinfo[[5]]][[1]]]+0.05width},{(Min[Transpose[Plotinfo[[5]]][[2]]]-0.03*height),(Max[Transpose[Plotinfo[[5]]][[2]]]+0.05*height)}},AspectRatio->Full,AxesLabel->{Plotinfo[[1,2]],"μ"<>ToString[StringDrop[Plotinfo[[4]],-2]]},PlotLabel->ToString[Plotinfo[[1,1]]]<>" Relative RMS ="<>ToString[N[Round[flagval*1000]/1000]],Epilog->Table[{Transpose[Plotinfo[[2]]][[2,i]],PointSize[0.01],Point[Plotinfo[[5,i]]],Text[Plotinfo[[2,i,1]],Plotinfo[[5,i]]+{0.0,0.03*height}]},{i,1,Length[Plotinfo[[2]]]}]]
];
MuAllele2FuncPlot[Plotinfo_,func_,flagval_]:=Block[{
height=(Max[Transpose[Plotinfo[[6]]][[2]]]-Min[Transpose[Plotinfo[[6]]][[2]]]),
width=(Max[Transpose[Plotinfo[[6]]][[1]]]-Min[Transpose[Plotinfo[[6]]][[1]]])},
Plot[func[var],{var,Min[Transpose[Plotinfo[[6]]][[1]]],Max[Transpose[Plotinfo[[6]]][[1]]]},
PlotStyle->AbsoluteThickness[Min[{(Maxflagvalue/flagval),10}]],PlotRange->{{Min[Transpose[Plotinfo[[6]]][[1]]]-0.05width,Max[Transpose[Plotinfo[[6]]][[1]]]+0.05width},{(Min[Transpose[Plotinfo[[6]]][[2]]]-0.03*height),(Max[Transpose[Plotinfo[[6]]][[2]]]+0.05*height)}},AspectRatio->Full,AxesLabel->{Plotinfo[[1,2]],"μ"<>ToString[StringDrop[Plotinfo[[4]],2]]},PlotLabel->ToString[Plotinfo[[1,1]]]<>" Relative RMS ="<>ToString[N[Round[flagval*1000]/1000]],Epilog->Table[{Transpose[Plotinfo[[2]]][[2,i]],PointSize[0.01],Point[Plotinfo[[6,i]]],Text[Plotinfo[[2,i,1]],Plotinfo[[6,i]]+{0.0,0.03*height}]},{i,1,Length[Plotinfo[[2]]]}]]
];

Maxflagvalue=0.1;  MinPops=7;
For[j=1,j≤Length[SNPsFileNames],j++,
For[k=1,k≤Length[EnvParDA[[1]]]-1,k++,
flagvalue1=MyFun[1][SNPAndAllelicMusvsParam[k,j][[3]],musfunc[(Cos[b]^2*(Tanh[a*(x/Sort[Transpose[SNPAndAllelicMusvsParam[k,j][[3]]][[1]]][[-1]])+x0]+1)/2+Sin[c]^2)/2,2*Tp,Mu],{x0,b,c,Tp,Mu,a},x];
flagvalue4=MyFun[4][SNPAndAllelicMusvsParam[k,j][[3]],a*x+b,{b,a},x];
flagvalue7=MyFun[7][SNPAndAllelicMusvsParam[k,j][[3]],a*(x^2)+b*x+c,{b,a,c},x];
flagvalue2=MyFun[2][SNPAndAllelicMusvsParam[k,j][[5]],muafunc[(Cos[b]^2*(Tanh[a*(x/Sort[Transpose[SNPAndAllelicMusvsParam[k,j][[5]]][[1]]][[-1]])+x0]+1)/2+Sin[c]^2)/2,2*Tp,Mu],{x0,b,c,Tp,Mu,a},x];
flagvalue5=MyFun[5][SNPAndAllelicMusvsParam[k,j][[5]],a*x+b,{b,a},x];
flagvalue8=MyFun[8][SNPAndAllelicMusvsParam[k,j][[5]],a*(x^2)+b*x+c,{b,a,c},x];
flagvalue3=MyFun[3][SNPAndAllelicMusvsParam[k,j][[6]],muafunc[(Cos[b]^2*(Tanh[a*(x/Sort[Transpose[SNPAndAllelicMusvsParam[k,j][[6]]][[1]]][[-1]])+x0]+1)/2+Sin[c]^2)/2,2*Tp,Mu],{x0,b,c,Tp,Mu,a},x];
flagvalue6=MyFun[6][SNPAndAllelicMusvsParam[k,j][[6]],a*x+b,{b,a},x];
flagvalue9=MyFun[9][SNPAndAllelicMusvsParam[k,j][[6]],a*(x^2)+b*x+c,{b,a,c},x];
AlleleOneValues=Sort[{flagvalue2,flagvalue5,flagvalue8}];
AlleleTwoValues=Sort[{flagvalue3,flagvalue6,flagvalue9}];
AllValuesSorted=Sort[{AlleleOneValues[[1]],AlleleTwoValues[[1]]}];
If[(AllValuesSorted[[2]]<=Maxflagvalue)&&(Length[SNPAndAllelicMusvsParam[k,j][[3]]]≥MinPops),
flagplots[k,j]={
SNPfuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[1],flagvalue1],
SNPfuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[4],flagvalue4],
SNPfuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[7],flagvalue7],
MuAllele1FuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[2],flagvalue2],
MuAllele1FuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[5],flagvalue5],
MuAllele1FuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[8],flagvalue8],
MuAllele2FuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[3],flagvalue3],
MuAllele2FuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[6],flagvalue6],
MuAllele2FuncPlot[SNPAndAllelicMusvsParam[k,j],Ftfunc[9],flagvalue9]
},##&[]
]
//Print
];
];