Clear["Global`*"];
Clear[model, data, ρeff, Dexp,temp, T, visc, MW, MWs, α];

nsolvents=5;

(*preparation for the presentation of data*)
imagesize={500,500};

square=Graphics[{EdgeForm[Thick],White,Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
squareR=Graphics[{Red, Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
squareR2=Graphics[{EdgeForm[Thin],LightRed,Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
squareO=Graphics[{Orange, Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
squareB=Graphics[{Black, Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
squareBl=Graphics[{Blue, Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
diamond=Graphics[{EdgeForm[Thick],White,Polygon[{{-.01,.0},{.0,.01},{.01,.0},{.0,-.01}}]}];
diamondBl=Graphics[{Black,Polygon[{{-.01,.0},{.0,.01},{.01,.0},{.0,-.01}}]}];
diamondP=Graphics[{Purple,Polygon[{{-.01,.0},{.0,.01},{.01,.0},{.0,-.01}}]}];
diamondB=Graphics[{Blue,Polygon[{{-.01,.0},{.0,.01},{.01,.0},{.0,-.01}}]}];
circle=Graphics[{EdgeForm[{Thick,Black}],Thickness[0.33],Circle[{0,0},1]}];
circleB=Graphics[{Black,Disk[{0,0},1]}];
triangle=Graphics[{EdgeForm[Thick],White,Polygon[{{0,1},{Sin[2 Pi/3],Cos[2 Pi/3]},{Sin[4 Pi/3],Cos[4 Pi/3]}}]}];
triangleB=Graphics[{Black,Polygon[{{0,1},{Sin[2 Pi/3],Cos[2 Pi/3]},{Sin[4 Pi/3],Cos[4 Pi/3]}}]}];
solidsquare=Graphics[Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]];
soliddiamond=Graphics[Polygon[{{-.01,.0},{.0,.01},{.01,.0},{.0,-.01}}]];

plotmarkersD={diamondB,0.025};
plotmarkersDRE={squareR,0.02};
plotmarkersD2={squareR2,0.02};
plotmarkersDRE2={squareB,0.02};
plotmarkersDvsD={{diamondB,0.025}, {squareR,0.02}};
plotmarkersCDCl3={diamond, 0.025};
plotmarkersD2O={square, 0.018};
plotmarkersTOL={circle, 0.02};
plotmarkersDMSO={triangle, 0.02};
plotmarkersPURE={squareB,0.02};
plotmarkersBENZ={triangleB, 0.02};
plotmarkersTHF={circleB, 0.02};
plotmarkersCYC={diamondBl, 0.02};
plotmarkersC298D2O={squareO,0.018};
plotmarkersC303D2O={squareR,0.018};
plotmarkersC298CDCl3={diamondB, 0.025};
plotmarkersC303CDCl3={diamondP, 0.025};
plotmarkersC={{squareO,0.018}, {squareR,0.018}, {diamondB, 0.025}, {diamondP, 0.025}};
plotmarkersV={{diamond,.023},{square,.023},{triangle,.03},{solidsquare,.022},{soliddiamond,.03}};

drange={{0,50},{0,50}};
deldrange={{0,1500},{-70,70}};
deldrange2={{0,50},{-70,70}};
dvrange={{0,25},{0,25}};
dvranges={{0,2.5 10^1},{0,2.5 10^1}};

slope45=Plot[Tooltip[x,"theoretical line"],{x,0,50},AspectRatio->1,ImageSize->imagesize,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.003],LabelStyle->18];
line25mw=Plot[Tooltip[25,"theoretical line"],{x,0,1500},AspectRatio->1,ImageSize->imagesize,PlotStyle->{Black,Dashed,Thickness[0.003]},AxesStyle->Thickness[0.003],LabelStyle->18];
linem25mw=Plot[Tooltip[-25,"theoretical line"],{x,0,1500},AspectRatio->1,ImageSize->imagesize,PlotStyle->{Black,Dashed,Thickness[0.003]},AxesStyle->Thickness[0.003],LabelStyle->18];
line25d=Plot[Tooltip[25,"theoretical line"],{x,0,50},AspectRatio->1,ImageSize->imagesize,PlotStyle->{Black,Dashed,Thickness[0.003]},AxesStyle->Thickness[0.003],LabelStyle->18];
linem25d=Plot[Tooltip[-25,"theoretical line"],{x,0,50},AspectRatio->1,ImageSize->imagesize,PlotStyle->{Black,Dashed,Thickness[0.003]},AxesStyle->Thickness[0.003],LabelStyle->18];


(* import data from Rob's Desktop*)
path="C:\\Users\\evansr2\\";
SetDirectory[path<>"\\Desktop"]; 
(* import data from GAM HDD 
SetDirectory["/Users/gam/Documents/WP/Papers/RE GDP D(MW) 2/Anal Chem revision"]*)
FileNames[];
Data=Import["Improving the Quantitative Interpretation of diffusion-ordered NMR spectra - all data.xlsx"];
Data=Partition[Flatten[Data], 15];

DataT={Data[[2]]};
For[i=3, i≤667,i++, DataT=AppendTo[DataT, Data[[ i]]]];

Dimensions[DataT];
entries=Dimensions[DataT, 1] ;

SeedRandom[];
loopsize=5;
For [randomloopcount=1,randomloopcount≤ loopsize,randomloopcount++,
DataRANDOM=RandomSample[DataT, 666]; 
Dimensions[DataRANDOM];

(*produce training set*)

(*original Dexp*)
Dexp={DataRANDOM[[1, 12]]};
For[i=2, i≤333,i++, Dexp=AppendTo[Dexp, DataRANDOM[[ i, 12]]]];
Dexp;

(*original T data*)
T={DataRANDOM[[1, 4]]};
For[i=2, i≤333,i++, T=AppendTo[T, DataRANDOM[[ i, 4]]]];
T;

(*original visc data*)
visc={DataRANDOM[[1, 8]]};
For[i=2, i≤333,i++, visc=AppendTo[visc, DataRANDOM[[ i, 8]]]];
visc;

(*original MW data*)
MW={DataRANDOM[[1, 6]]};
For[i=2, i≤333,i++, MW=AppendTo[MW, DataRANDOM[[ i, 6]]]];
MW;

(*original MWs data*)
MWs={DataRANDOM[[1, 7]]};
For[i=2, i≤333,i++, MWs=AppendTo[MWs, DataRANDOM[[ i, 7]]]];
MWs;

α={(MWs[[1]]/MW[[1]])^(1/3)};
For[i=2, i≤333,i++, α=AppendTo[α, (MWs[[i]]/MW[[i]])^(1/3)]];

data={{MW[[1]], MWs[[1]], visc[[1]], T[[1]], Dexp[[1]]}};
For[i=2, i≤333,i++,data=AppendTo[data, { MW[[i]], MWs[[i]], visc[[i]],T[[i]], Dexp[[i]]}]];
data;
nd=Length[data];

Clear[model];
α=(mws/mw)^(1/3);
model=10^10(kB temp((1.5α)+(1/(1+α))))/(6 Pi η (0.003mw/(4000 Pi ρeff nA))^(1/3));
(* Print["original model"]
model *)

kB=1.380662 10^-23;
nA=6.022045 10^23;

Print["minimising (Dexp-Dcalc)/Dcalc) for the training data set of 333 randomly chosen measurements gives"];

model2=dexp/model;
data2=Transpose[Join[Transpose[data],{Table[1,{j,Length[data]}]}]];

nlm=NonlinearModelFit[data2,model2,ρeff,{mw,mws,η,temp,dexp}]; (* minimise |(Dexp-Dcalc)/Dcalc| *)
(* nlm["ParameterConfidenceIntervals"] for a single variable, CI is ± 1.96 standard error *)
ρeffoptorig=ρeff/.nlm["BestFitParameters"]; ρefferr=nlm["ParameterErrors"][[1]]; 
Print["ρeff = ",ρeffoptorig," ± 1.96 * ",ρefferr];


(*create training set of original data*)
Dcalctraining = {{model/.{mw-> data[[1,1]] ,mws-> data[[1,2]],η->  data[[1,3]], temp->  data[[1,4]],ρeff->ρeffoptorig},data[[1,5]]}};

For[i=2,i≤333,i++, AppendTo[Dcalctraining ,{model/.{mw-> data[[i,1]] ,mws-> data[[i,2]],η->  data[[i,3]], temp->  data[[i,4]], ρeff->ρeffoptorig},data[[i,5]]}];];

dratios1=Table[Dcalctraining [[i,2]]/Dcalctraining [[i,1]],{i,1,Length[Dcalctraining ]}];

(*produce test set*)

(*original Dexp*)
testDexp={DataRANDOM[[334, 12]]};
For[i=335, i≤666,i++, testDexp=AppendTo[testDexp, DataRANDOM[[ i, 12]]]];
testDexp;

(*original T data*)
testT={DataRANDOM[[334, 4]]};
For[i=335, i≤666,i++, testT=AppendTo[testT, DataRANDOM[[ i, 4]]]];
testT;

(*original visc data*)
testvisc={DataRANDOM[[334, 8]]};
For[i=335, i≤666,i++, testvisc=AppendTo[testvisc, DataRANDOM[[ i, 8]]]];
testvisc;

(*original MW data*)
testMW={DataRANDOM[[334, 6]]};
For[i=335, i≤666,i++, testMW=AppendTo[testMW, DataRANDOM[[ i, 6]]]];
testMW;

(*original MWs data*)
testMWs={DataRANDOM[[334, 7]]};
For[i=335, i≤666,i++, testMWs=AppendTo[testMWs, DataRANDOM[[ i, 7]]]];
testMWs;

testα={(testMWs[[1]]/testMW[[1]])^(1/3)};
For[i=2, i≤333,i++, testα=AppendTo[testα, (testMWs[[i]]/testMW[[i]])^(1/3)]];

testdata={{testMW[[1]], testMWs[[1]], testvisc[[1]], testT[[1]], testDexp[[1]]}};
For[i=2, i≤333,i++,testdata=AppendTo[testdata, { testMW[[i]], testMWs[[i]], testvisc[[i]],testT[[i]], testDexp[[i]]}]];
testdata;
testnd=Length[testdata];

Clear[mw, mws, η, temp,ρeff];

(*create test set of original data*)
Dcalctest= {{model/.{mw->testdata[[1,1]]  ,mws->testdata[[1,2]] ,η->  testdata[[1,3]] , temp->  testdata[[1,4]],ρeff->ρeffoptorig},testdata[[1,5]]}};
For[i=2,i≤333,i++, AppendTo[Dcalctest,{model/.{mw->testdata[[i,1]]  ,mws->testdata[[i,2]] ,η-> testdata[[i,3]] ,temp->  testdata[[i,4]], ρeff->ρeffoptorig},testdata[[i,5]]}];];

Print["minimising (Dexp-Dcalc)/Dcalc) for the remaining test data set of 333 randomly chosen measurements gives"];

model2=dexp/model;
testdata2=Transpose[Join[Transpose[testdata],{Table[1,{j,Length[testdata]}]}]];

nlm=NonlinearModelFit[testdata2,model2,ρeff,{mw,mws,η,temp,dexp}]; (* minimise |(Dexp-Dcalc)/Dcalc| *)
(* nlm["ParameterConfidenceIntervals"] for a single variable, CI is ± 1.96 standard error *)
ρeffopttest=ρeff/.nlm["BestFitParameters"]; ρefferr=nlm["ParameterErrors"][[1]]; 
Print["ρeff = ",ρeffopttest," ± 1.96 * ",ρefferr];

dratios=Table[Dcalctest[[i,2]]/Dcalctest[[i,1]],{i,1,Length[Dcalctest]}];

Print["RMS percentage deviation from ratio fit model using density optimised from training data set, applied to training data set ",100 (Sum[(dratios1[[i]]-1)^2,{i,1,Length[Dcalctraining ]}]/Length[Dcalctraining ])^(1/2)];Print["RMS percentage deviation from ratio fit model using density optimised from training data set, applied to test data set ",100 (Sum[(dratios[[i]]-1)^2,{i,1,Length[Dcalctest]}]/Length[Dcalctest])^(1/2)];

(*create test set of original data*)
Dcalctest2= {{model/.{mw->testdata[[1,1]]  ,mws->testdata[[1,2]] ,η->  testdata[[1,3]] , temp->  testdata[[1,4]],ρeff->ρeffopttest},testdata[[1,5]]}};
For[i=2,i≤333,i++, AppendTo[Dcalctest2,{model/.{mw->testdata[[i,1]]  ,mws->testdata[[i,2]] ,η-> testdata[[i,3]] ,temp->  testdata[[i,4]], ρeff->ρeffopttest},testdata[[i,5]]}];];

dratios2=Table[Dcalctest2[[i,2]]/Dcalctest2[[i,1]],{i,1,Length[Dcalctest]}];

Print["RMS percentage deviation from ratio fit model using density optimised from test data set, applied to test data set ",100 (Sum[(dratios2[[i]]-1)^2,{i,1,Length[Dcalctest2]}]/Length[Dcalctest2])^(1/2)]
];
Print[];
Print["Dexp vs Dcalc, Dcalc estimated using density produced from training set and applied to test data set of remaining 333 measurements"];

DvsD=ListPlot[Dcalctest, PlotRange->drange,AspectRatio->0.88,ImageSize->imagesize, PlotMarkers->plotmarkersD,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.01],Frame->{{True,True},{True,True}},FrameStyle->{{Thick, Thick},{Thick, Thick}},FrameLabel->{Style["\!\(\*TagBox[SubscriptBox[
  StyleBox[
   StyleBox[\"D\",\nStripOnInput->False,\nFontSlant->Italic],\nStripOnInput->False,\nFontSlant->Italic], \"est\"],
 DisplayForm]\) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False],Style["\!\(\*TagBox[SubscriptBox[
  StyleBox[
   StyleBox[\"D\",\nStripOnInput->False,\nFontSlant->Italic],\nStripOnInput->False,\nFontSlant->Italic], \"exp\"],
 DisplayForm]\) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False]},LabelStyle->{24,Black, FontFamily->"Helvetica"}];

Show[DvsD, slope45]