(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
BeginPackage["DataAnalysis`CompositionalAnalysis`",{
"SciDraw`",
"DataAnalysis`"
}];


(* ::Input::Initialization:: *)
Unprotect[Evaluate[$Context<>"*"]];


(* ::Input::Initialization:: *)
CalibrationPlot;
CalibrateCompositionalData;
RemoveNonDetectableValues;

FilterOutliers;

SSBSubstractBlanks;
SSBCalculateDeltas;
SSBCalculateMeanRatios;
SSBCalculateCalibrationFit;


(* ::Input::Initialization:: *)
IncludeQuadraticTerm;
DriftStandard;
BlankPattern;


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SetAttributes[ErrorMessage,HoldFirst];
ErrorMessage[Args___]:=Module[
{},
Message[Args];
Abort[];
];


(* ::Input::Initialization:: *)
SetAttributes[CatchError,HoldAll];
CatchError[Body_]:=(Catch[Body]===$ErrorFlag);


(* ::Input::Initialization:: *)
Options[CalibrationPlot]={IncludeConstantBasis->False,IncludeQuadraticTerm->False,BlankPattern->Null};
BlankPattern::usage="BlankPattern is an option to specify a string pattern that matches the row keys (IDs) of all blank measurements.";

CalibrationPlot::nostandards="No standards were found for this element in the provided standard data.";
CalibrationPlot::nostandarddata="No data was found for these standards and element in the provided input data.";
CalibrationPlot::noisotope="No element/isotope with the provided ID was measured.";

CalibrationPlot::usage="CalibrationPlot[datatable,standardData,column,{\!\(\*SubscriptBox[\(standardPatttern\), \(1\)]\),\!\(\*SubscriptBox[\(standardPatttern\), \(2\)]\),...}] plots the measurements in a specific column (e.g. compound) of the DataTable for the standards matched by each \!\(\*SubscriptBox[\(standardPatttern\), \(i\)]\) vs. the corresponding concentrations specified by the standards.
CalibrationPlot[datatable,standardData,{\!\(\*SubscriptBox[\(column\), \(i\)]\),\!\(\*SubscriptBox[\(column\), \(2\)]\),...},{\!\(\*SubscriptBox[\(standardPatttern\), \(1\)]\),\!\(\*SubscriptBox[\(standardPatttern\), \(2\)]\),...}] plots the measurements in a series of columns (e.g. compounds).";
CalibrationPlot[intData_DataTable,standardConc_DataTable,column_String,stdPatterns_List,opts:OptionsPattern[]]:=Module[
{FigDat,YMax,XMax,FitF},
(
FigDat=CalibrationPlotData[intData,standardConc,column,#,opts]&/@stdPatterns;
FitF=FigDat[[All,2]];
FigDat=FigDat[[All,1]];
FigDat=DeleteMissing[#,1,Infinity]&/@FigDat;

XMax=Max[(Max@#[All,"Concentration",DVMeans])&/@FigDat];
XMax=Ceiling[XMax,10^((Log10[XMax]//Floor))];

YMax=Max[(Max@#[All,"Measurement",DVMeans])&/@FigDat];
YMax=Ceiling[YMax,10^((Log10[YMax]//Floor))];

Show[
DTPlot[
FigDat,
PlotRange->{{0,XMax},{0,YMax}},
PlotLegends->stdPatterns,
PlotLabel->Style[column,16,Black],
ImageSize->500,
PlotRangePadding->Scaled[0.05]
],
MapIndexed[Plot[If[TrueQ[OptionValue[IncludeQuadraticTerm]],
If[TrueQ[OptionValue[IncludeConstantBasis]],
a*x^2+b*x+c,
a*x^2+b*x
],
If[TrueQ[OptionValue[IncludeConstantBasis]],
b*x+c,
b*x
]
]/.#["BestFitParameters"],{x,0,XMax},PlotStyle->Black]&,FitF]
]

(*DataAnalysis`Fig[
Join[
MapIndexed[Plot[If[TrueQ[OptionValue[IncludeQuadraticTerm]],
If[TrueQ[OptionValue[IncludeConstantBasis]],
a*x^2+b*x+c,
a*x^2+b*x
],
If[TrueQ[OptionValue[IncludeConstantBasis]],
b*x+c,
b*x
]
]/.#["BestFitParameters"],{x,0,XMax},PlotStyle\[Rule]Black]&,FitF],
FigDat
]
]*)
)
];

CalibrationPlot[intData_DataTable,standardConc_DataTable,columns_List,stdPatterns_List,opts:OptionsPattern[]]:=Module[
{PlotDat,FigDat,YMax,XMax,FitF},
(
Column[(CalibrationPlot[intData,standardConc,#,stdPatterns,opts]&/@columns),Center,5]
)
]


(* ::Input::Initialization:: *)
Options[CalibrationPlotData]={IncludeConstantBasis->False,IncludeQuadraticTerm->False,BlankPattern->Null};
CalibrationPlotData[intData_DataTable,standardConc_DataTable,column_String,stdPattern_,opts:OptionsPattern[]]:=Module[
{Element,StandardConc,StandardInt,Standards,FitDat,FitDatE,FitF,IDat,OutDat,StandardsFound,FD},
(
Off[InterpolatingFunction::dmval];
Off[FittedModel::varzero];
Off[Interpolation::inhr];

Element=StringCases[column,LetterCharacter..]//First;
StandardConc=standardConc[All,Element~~___];

If[FailureQ[StandardConc]||(StandardConc[All,Length][Max]==0),ErrorMessage[CalibrationPlot::nostandards]];
(*StandardConc=ReplaceRight[StandardConc,(Select[Transpose[#],DVMeans[#[[2]]>0]&]//Transpose)&];*)
Standards=Normal@Keys@StandardConc;

If[Length[intData[stdPattern]]==0,ErrorMessage[CalibrationPlot::nostandarddata]];
StandardInt=Flatten@Map[Normal@Values@intData[(___~~#~~___),column][KeyFind[stdPattern]]&,Standards];
If[MatchQ[StandardInt,{Null..}],ErrorMessage[CalibrationPlot::noisotope]];

StandardsFound=DeleteCases[Map[If[Length[intData[(___~~#~~___),column][KeyFind[stdPattern]]]==1,#,Null]&,Standards],Null];

IDat=CalibrationInterpolations[intData,StandardsFound,column,opts];

FitDat={
DVMeans@Flatten@Values@Values@Normal@StandardConc[StandardsFound],
(StandardInt//DVMeans)-(IDat[[3]][#]&/@(Values@Normal@intData[stdPattern,AbsoluteTime[#"Time"]&]))
}//Transpose;

FitDatE={
DVSEs@Flatten@Values@Values@Normal@StandardConc[StandardsFound],
(StandardInt//DVSEs)+(IDat[[4]][#]&/@(Values@Normal@intData[stdPattern,AbsoluteTime[#"Time"]&]))
}//Transpose;

FD=MapThread[{#1,#2}&,{FitDat,FitDatE}];
FD=DeleteMissing[FD,1,Infinity];

FitF=If[TrueQ[OptionValue[IncludeQuadraticTerm]],
If[TrueQ[OptionValue[IncludeConstantBasis]],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],a*x^2+b*x+c,{a,b,c},{x},5],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],a*x^2+b*x,{a,b},{x},5]
],
If[TrueQ[OptionValue[IncludeConstantBasis]],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],b*x+c,{b,c},{x},5],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],b*x,{b},{x},5]
]
];

OutDat=DataTable[Prepend[MapThread[{#3,#1,DVE@(#2-NDV[IDat[[3]][#4],IDat[[4]][#4]])}&,{Flatten@Values@Values@Normal@StandardConc[StandardsFound],StandardInt,StandardsFound,(Values@Normal@intData[stdPattern,AbsoluteTime[#"Time"]&])}],{"ID","Concentration","Measurement"}]];

On[InterpolatingFunction::dmval];
On[FittedModel::varzero];
On[Interpolation::inhr];

{OutDat,FitF}
)
]


(* ::Input::Initialization:: *)
Options[CalibrateCompositionalData]={IncludeConstantBasis->False,IncludeQuadraticTerm->False,BlankPattern->Null,DriftStandard->Null};
DriftStandard::usage="DriftStandard is an option to specify a string pattern that matches the row keys (IDs) of drift standard measurements.";

CalibrateCompositionalData::nostandards="No standards were found for this element in the provided standard data.";
CalibrateCompositionalData::noisotope="No element/isotope with the provided ID was measured.";

CalibrateCompositionalData::usage="CalibrateCompositionalData[datatable,standardData,{\!\(\*SubscriptBox[\(column\), \(i\)]\),\!\(\*SubscriptBox[\(column\), \(2\)]\),...}] calibrates the measurements for the compounds in a set of columns in a DataTable using the standardData.";
CalibrateCompositionalData[intData_DataTable,standardConc_DataTable,columns_List,opts:OptionsPattern[]]:=Module[
{IntDt,Element,StandardConc,Standards,StandardInt,StandardsFound,OutDt,IDat,BInterpols={},BInterpolsE={},Interpols={},InterpolsE={},DriftInterpols={},DriftInterpolsE={},StandardConcs={},ColName,ColNum,FullColNames={}},
(
IntDt=intData;

Off[InterpolatingFunction::dmval];
Off[FittedModel::varzero];
Off[Interpolation::inhr];

Scan[(
Element=StringCases[#,LetterCharacter..]//First;

StandardConc=standardConc[All,Element~~___];
If[FailureQ[StandardConc]||(StandardConc[All,Length][Max]==0),ErrorMessage[CalibrateCompositionalData::nostandards]];
Standards=Normal@Keys@StandardConc;

StandardInt=Flatten@Map[Normal@Values@intData[(___~~#~~___),column]&,Standards];
If[MatchQ[StandardInt,{Null..}],ErrorMessage[CalibrateCompositionalData::noisotope]];

StandardsFound=DeleteCases[Map[If[Length[intData[(___~~#~~___),column]]>=1,#,Null]&,Standards],Null];

IDat=CalibrationInterpolations[intData,StandardsFound,#,opts];

AppendTo[Interpols,IDat[[1]]];
AppendTo[InterpolsE,IDat[[2]]];
AppendTo[BInterpols,IDat[[3]]];
AppendTo[BInterpolsE,IDat[[4]]];
AppendTo[DriftInterpols,IDat[[5]]];
AppendTo[DriftInterpolsE,IDat[[6]]];
AppendTo[StandardConcs,StandardConc[StandardsFound]];
AppendTo[FullColNames,If[Length[StringCases[First@ColNames@StandardConc[StandardsFound]," / "~~__]]>0,#<>(StringCases[First@ColNames@StandardConc[StandardsFound]," / "~~__]//First),#]]
)&,columns];

OutDt=intData;
MapIndexed[
(
ColName=#1;
ColNum=#2[[1]];

QueryTo[OutDt,
All,(#~Join~Association[
{
(FullColNames[[ColNum]])->CalibrateDatum[#[ColName],ColName,StandardConcs[[ColNum]],#"Time",Interpols[[ColNum]],InterpolsE[[ColNum]],BInterpols[[ColNum]],BInterpolsE[[ColNum]],DriftInterpols[[ColNum]],DriftInterpolsE[[ColNum]],opts]
}
])&,
FailureAction->None
];
)&,columns
];
OutDt=OutDt/.{x_Complex->\[Placeholder]};

On[Interpolation::inhr];
On[InterpolatingFunction::dmval];
On[FittedModel::varzero];

OutDt[All,FullColNames]
)
];


(* ::Input::Initialization:: *)
Options[CalibrationInterpolations]={IncludeConstantBasis->False,IncludeQuadraticTerm->False,BlankPattern->Null,DriftStandard->Null};
CalibrationInterpolations[intData_,standardsFound_,column_,opts:OptionsPattern[]]:=Module[
{(*BInterpols,BInterpolsE,*)Interpols,InterpolsE,DriftReference,DriftInterpols=1&,DriftInterpolsE=0&,V},
(
(* Interpolations for standard intensities *)
Interpols=Map[(
V=Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][___~~#~~___,{AbsoluteTime[#"Time"],DVMeans@Slot[1][column]}&];
If[Length[V]>0,Interpolation[V,InterpolationOrder->1],Missing[]&]
)&,standardsFound];
InterpolsE=Map[
(
V=Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][___~~#~~___,{AbsoluteTime[#"Time"],DVSEs@Slot[1][column]}&];
If[Length[V]>0,Interpolation[V,InterpolationOrder->1],Missing[]&]
)&,standardsFound];

(* Interpolations for blank intensities *)
(*If[TrueQ[OptionValue[BlankPattern]\[Equal]Null],
BInterpols=(0&);
BInterpolsE=(0&);
,
BInterpols=Interpolation[Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][OptionValue[BlankPattern],{AbsoluteTime[#"Time"],DVMeans@Slot[1][column]}&],InterpolationOrder\[Rule]1];
BInterpolsE=Interpolation[Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][OptionValue[BlankPattern],{AbsoluteTime[#"Time"],DVSEs@Slot[1][column]}&],InterpolationOrder\[Rule]1];
];*)

If[StringQ[OptionValue[DriftStandard]],
V=Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][___~~(OptionValue[DriftStandard])~~___,{AbsoluteTime[#"Time"],DVMeans@Slot[1][column]}&];
DriftReference=If[Length[V]>0,Interpolation[V,InterpolationOrder->1],Missing[]&];

V=Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][___~~("DRIFT"|OptionValue[DriftStandard])~~___,{AbsoluteTime[#"Time"],(DVMeans@Slot[1][column])/DriftReference[AbsoluteTime[#"Time"]]}&,FailureAction->None];
DriftInterpols=If[Length[V]>0,Interpolation[V,InterpolationOrder->1],Missing[]&];

V=Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][___~~("DRIFT"|OptionValue[DriftStandard])~~___,{AbsoluteTime[#"Time"],(DVSEs@Slot[1][column])/DriftReference[AbsoluteTime[#"Time"]]}&,FailureAction->None];
DriftInterpolsE=If[Length[V]>0,Interpolation[V,InterpolationOrder->1],Missing[]&];
];

Join[{Interpols,InterpolsE},BlankInterpolations[intData,column,FilterRules[{opts}, Options[BlankInterpolations]]],{DriftInterpols,DriftInterpolsE}]
(*{Interpols,InterpolsE,BInterpols,BInterpolsE,DriftInterpols,DriftInterpolsE}*)
)
];


(* ::Input::Initialization:: *)
Options[BlankInterpolations]={BlankPattern->Null};
BlankInterpolations[intData_,column_,opts:OptionsPattern[]]:=Module[
{BInterpols,BInterpolsE},
(
If[TrueQ[OptionValue[BlankPattern]==Null],
BInterpols=(0&);
BInterpolsE=(0&);
,
BInterpols=Interpolation[Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][OptionValue[BlankPattern],{AbsoluteTime[#"Time"],DVMeans@Slot[1][column]}&],InterpolationOrder->1];
BInterpolsE=Interpolation[Values@Normal@intData[Select[NumericDVQ[Slot[1][column]]&]][OptionValue[BlankPattern],{AbsoluteTime[#"Time"],DVSEs@Slot[1][column]}&],InterpolationOrder->1];
];

{BInterpols,BInterpolsE}
)
]


(* ::Input::Initialization:: *)
Options[RemoveNonDetectableValues]={BlankPattern->Null};
RemoveNonDetectableValues::usage="RemoveNonDetectableValues[datatable,{\!\(\*SubscriptBox[\(column\), \(i\)]\),\!\(\*SubscriptBox[\(column\), \(2\)]\),...}] removes values in intData that are below the detection limit based on blank measurements specified by the BlankPattern option.";
RemoveNonDetectableValues[intData_DataTable,columns_List,opts:OptionsPattern[]]:=Module[
{IntDt,OutDt,BlankValues,ColName,ColNum},
(
IntDt=intData;

Off[InterpolatingFunction::dmval];
Off[Interpolation::inhr];

OutDt=intData;
MapIndexed[
(
ColName=#1;
ColNum=#2[[1]];

BlankValues=Values@Normal@intData[Select[NumericDVQ[Slot[1][ColName]]&]][OptionValue[BlankPattern],DVMeans@Slot[1][ColName]&];

QueryTo[OutDt,
KeyValueMap[<|#1->If[StringMatchQ[#1,OptionValue[BlankPattern]],#2,(#2~Join~Association[
{
ColName->TestDetectionLimit[#2[ColName],#2["Time"],Mean[BlankValues],StandardDeviation[BlankValues]]
}
])]|>&,#]&,
FailureAction->None
];

OutDt=DataTable[<|OutDt//Normal|>];
)&,columns
];

On[Interpolation::inhr];
On[InterpolatingFunction::dmval];

DataTable[<|OutDt//Normal|>]
)
]


(* ::Input::Initialization:: *)
TestDetectionLimit[datum_,time_,bMean_,bSTD_]:=Module[{},(
If[MissingQ[datum],Return[datum]];
If[!NumericQ[DVMeans[datum]],Return[datum]];

If[DVMeans[datum]<(bMean+3*bSTD),
Missing["BelowDL",(bMean+3*bSTD)],
(*If[DVMeans[datum]<(bInterpol[AbsoluteTime[time]]+10*bInterpolE[AbsoluteTime[time]]),
Missing["BelowQL",(bInterpol[AbsoluteTime[time]]+10*bInterpolE[AbsoluteTime[time]])],
datum
]*)
datum
]
)]


(* ::Input::Initialization:: *)
Options[CalibrateDatum]={IncludeConstantBasis->False,IncludeQuadraticTerm->False,BlankPattern->Null,DriftStandard->Null};
CalibrateDatum[datum_,column_,standardConc_,time_,interpol_,interpolE_,bInterpol_,bInterpolE_,driftInterpol_,driftInterpolE_,opts:OptionsPattern[]]:=Module[
{FitF,FitDat,FitDatE,WeightDat,Params,BlankEV,DriftEV,Output,FD,Val},
(
If[MatchQ[datum,Missing[Except["BelowDL"|"BelowQL"],___]],Return[datum]];

FitDat={
DVMeans@Flatten@Values@Values@Normal@standardConc,
#[AbsoluteTime[time]]&/@interpol-bInterpol[AbsoluteTime[time]]
}//Transpose;

FitDatE={
DVSEs@Flatten@Values@Values@Normal@standardConc,
#[AbsoluteTime[time]]&/@interpolE+bInterpolE[AbsoluteTime[time]]
}//Transpose;

FD=MapThread[{#1,#2}&,{FitDat,FitDatE}];
FD=DeleteMissing[FD,1,Infinity];

FitF=If[TrueQ[OptionValue[IncludeQuadraticTerm]],
If[TrueQ[OptionValue[IncludeConstantBasis]],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],a*x^2+b*x+c,{a,b,c},{x},5],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],a*x^2+b*x,{a,b},{x},5]
],
If[TrueQ[OptionValue[IncludeConstantBasis]],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],b*x+c,{b,c},{x},5],
EffectiveVarianceFit[FD[[All,1]],FD[[All,2]],b*x,{b},{x},5]
]
];

Params=If[TrueQ[OptionValue[IncludeQuadraticTerm]],
If[TrueQ[OptionValue[IncludeConstantBasis]],
If[Length[FD[[All,1]]]==3,
{
a->NDV[a/.FitF["BestFitParameters"],0],
b->NDV[b/.FitF["BestFitParameters"],0],
c->NDV[c/.FitF["BestFitParameters"],0]
},
{
a->NDV[a/.FitF["BestFitParameters"],FitF["ParameterErrors"][[1]]],
b->NDV[b/.FitF["BestFitParameters"],FitF["ParameterErrors"][[2]]],
c->NDV[c/.FitF["BestFitParameters"],FitF["ParameterErrors"][[3]]]
}
],
If[Length[FD[[All,1]]]==2,
{
a->NDV[a/.FitF["BestFitParameters"],0],
b->NDV[b/.FitF["BestFitParameters"],0]
},
{
a->NDV[a/.FitF["BestFitParameters"],FitF["ParameterErrors"][[1]]],
b->NDV[b/.FitF["BestFitParameters"],FitF["ParameterErrors"][[2]]]
}
]
],
If[TrueQ[OptionValue[IncludeConstantBasis]],
If[Length[FD[[All,1]]]==2,
{
b->NDV[b/.FitF["BestFitParameters"],0],
c->NDV[c/.FitF["BestFitParameters"],0]
},
{
b->NDV[b/.FitF["BestFitParameters"],FitF["ParameterErrors"][[1]]],
c->NDV[c/.FitF["BestFitParameters"],FitF["ParameterErrors"][[2]]]
}
],
If[Length[FD[[All,1]]]==1,
{
b->NDV[b/.FitF["BestFitParameters"],0]
},
{
b->NDV[b/.FitF["BestFitParameters"],FitF["ParameterErrors"][[1]]]
}
]
]
];

BlankEV=If[bInterpolE[AbsoluteTime[time]]>0,NDV[bInterpol[AbsoluteTime[time]],bInterpolE[AbsoluteTime[time]]],bInterpol[AbsoluteTime[time]]];
DriftEV=If[NumericQ[driftInterpolE[AbsoluteTime[time]]],If[driftInterpolE[AbsoluteTime[time]]>0,NDV[driftInterpol[AbsoluteTime[time]],driftInterpolE[AbsoluteTime[time]]],driftInterpol[AbsoluteTime[time]]],1];

Val=datum;
If[MissingQ[datum],(Val=datum[[2]])];

Output=If[TrueQ[OptionValue[IncludeQuadraticTerm]],
If[TrueQ[OptionValue[IncludeConstantBasis]],
(
(-b+Sqrt[b^2+4*a*((Val/DriftEV-BlankEV))-4*a*c])/(2*a)
),
(-b+Sqrt[b^2+4*a*((Val/DriftEV-BlankEV))])/(2*a)
],
If[TrueQ[OptionValue[IncludeConstantBasis]],
(-c+((Val/DriftEV-BlankEV)))/b,
((Val/DriftEV-BlankEV))/b
]
]/.Params;
Output=DVEvaluate[Output];

If[MissingQ[datum],Missing[datum[[1]],Output],Output]
)
];


(* ::Input::Initialization:: *)
FilterOutliers::usage="FilterOutliers[data,{\!\(\*SubscriptBox[\(column\), \(1\)]\),\!\(\*SubscriptBox[\(column\), \(2\)]\),...},i] filters outliers in the specified columns (or keys) in data (DataTable or Association) based on a difference between the values and their mean which is larger than i standard deviations.
FilterOutliers[data,{\!\(\*SubscriptBox[\(column\), \(1\)]\),\!\(\*SubscriptBox[\(column\), \(2\)]\),...}] filters outliers based on differences larger than 2 standard devations.
FilterOutliers[data] filters outliers in all columns based on differences larger than 2 standard devations.";
FilterOutliers[data_DataTable,columns_:All,SDMultiple_:2]:=Module[
{EVs,Cols},
(
Cols=If[columns===All,ColNames@data,If[ListQ[columns],columns,List[columns]]];
EVs=data[All,Cols][MeanAndStandardDeviation];

data[Select[And@@(Module[{RD=#},(
(Abs[RD[#]-EVs[#][[1]]]<=SDMultiple*EVs[#][[2]])&/@Cols
)])&]]
)
];

FilterOutliers[data_Association,columns_:All,SDMultiple_:2]:=Normal@FilterOutliers[DataTable[data],columns];


(* ::Input::Initialization:: *)
SSBSubstractBlanks::usage="SSBSubstractBlanks[datatable] substracts the blank measurements from the standard and sample measurements in data which was obtained using the standard-sample-bracketing technique.";
SSBSubstractBlanks[data_DataTable]:=data[Partition[Normal@#,3,2]&][Association@#&,Module[
{Blanc},
(
Blanc=Merge[{Merge[Values@Values @#[[1]],Mean],Merge[Values@Values @#[[3]],Mean]},Mean];
Keys[#[[2]]]->(#-(Blanc)&/@Values@#[[2]])
)
]&];


(* ::Input::Initialization:: *)
RoughnessFun[deg_,knots_]:=RoughnessFun[deg,knots]=(
Integrate[D[Total@MapIndexed[ParamVar[#2[[1]]]*#1&,Array[XVar^#&,deg+1,0]~Join~Table[BSplineBasis[{deg,knots},i,XVar],{i,0,Length[knots]-deg-2}]],{XVar,2}]^2,{XVar,Min[knots],Max[knots]}]
);


(* ::Input::Initialization:: *)
RoughnessFun[deg_,knots_,params:{_?NumericQ...}]:=(
NIntegrate[D[Total@MapIndexed[params[[#2[[1]]]]*#1&,Array[XVar^#&,deg+1,0]~Join~Table[BSplineBasis[{deg,knots},i,XVar],{i,0,Length[knots]-deg-2}]],{XVar,2}]^2,{XVar,Min[knots],Max[knots]}]
);


(* ::Input::Initialization:: *)
CubicSplSmooth[data_,lambda_]:=Module[{M,Knots,X,Dsq,a},
M=Length@data;
Knots=Flatten@{1,1,1,Range@M,M,M,M};
X=Table[Evaluate@N@BSplineBasis[{3,Knots},n,t],{t,1,M},{n,0,M+1}];

Dsq=Differences[X,2];
a=LinearSolve[Transpose[X].X+lambda*Transpose[Dsq].Dsq,Transpose[X].data,Method->"Multifrontal"];
Return[X.a]
];


(* ::Input::Initialization:: *)
SmoothSplineFit[data_,errors_,deg_,lambda_]:=Module[
{Basis,BasisFun,SquaredErrors,ParamList,FitRes},
(
Basis=Array[XVar^#&,deg+1,0]~Join~Table[BSplineBasis[{deg,data[[All,1]]},i,XVar],{i,0,Length[data]-deg-2}];
BasisFun=Total@MapIndexed[ParamVar[#2[[1]]]*#1&,Basis];
ParamList=Array[ParamVar,Length[Basis]];
(*SquaredErrors=Simplify[Total[(1/#[[2]]^2)*(#[[1,2]]-(BasisFun/.XVar\[Rule]#[[1,1]]))^2&/@MapThread[{#1,#2}&,{data,errors}]]];
ParamList=Array[ParamVar,Length[Basis]];

FitRes=NMinimize[
(p)*SquaredErrors+(1-p)*RoughnessFun[deg,data[[All,1]]],
ParamList,
PrecisionGoal\[Rule]1,
AccuracyGoal\[Rule]1
];*)

FitRes=LinearModelFit[CubicSplSmooth[data,lambda],Basis,XVar];

BasisFun/.MapThread[#1->#2&,{ParamList,FitRes["BestFitParameters"]}]
)
]


(* ::Input::Initialization:: *)
SmoothSplines[data_,errors_,deg_,knots:{_?NumericQ..},rescaling_]:=Module[
{Basis,ExactFit,ParamList,BasisFun,SquaredErrors,StartParameters,SmoothFit,SquaredErrorFun,JacknifeScatter,ODFit},
(
(*Basis=Array[XVar^#&,deg+1,0]~Join~Table[BSplineBasis[{deg,knots},i,XVar],{i,0,Length[knots]-deg-2}];
ExactFit=LinearModelFit[data,Basis,XVar];

ParamList=Array[ParamVar,Length[Basis]];
BasisFun=Total@MapIndexed[ParamVar[#2[[1]]]*#1&,Basis];*)

(*SquaredErrors=Simplify[Total[(1/#[[2]]^2)*(#[[1,2]]-(BasisFun/.XVar\[Rule]#[[1,1]]))^2&/@MapThread[{#1,#2}&,{data,errors}]]];
*)
(*StartParameters=MapThread[{#1,#2}&,{ParamList,ExactFit["BestFitParameters"]}];*)

(*SquaredErrorFun[p_?NumericQ]:=SquaredErrorFun[p]=Module[{FitRes},
(
FitRes=NMinimize[
(p)*SquaredErrors+(1-p)*RoughnessFun[deg,knots],
ParamList,
PrecisionGoal\[Rule]2,
AccuracyGoal\[Rule]2
];

(SquaredErrors/.FitRes[[2]])
)
];*)

JacknifeScatter[lambda_?NumericQ]:=JacknifeScatter[lambda]=Total[
(
((SmoothSplineFit[Delete[data,#],Delete[errors,#],deg,lambda]/.XVar->data[[#,1]])-data[[#,2]])^2
)&/@Range[1,Length[data]]
];

SmoothFit=NMinimize[
{
JacknifeScatter[lambda],
lambda>0
},
{lambda},
Method->{"NelderMead","PostProcess"->False},
AccuracyGoal->3,
PrecisionGoal->3(*,
EvaluationMonitor\[RuleDelayed]Print[{lambda,JacknifeScatter[lambda]}]*)
];

SmoothFit=SmoothSplineFit[data,errors,3,lambda/.SmoothFit[[2]]];

ODFit=NMinimize[
{
(Mean[(1/(#[[2]]^2+od^2))*(#[[1,2]]-(SmoothFit/.XVar->#[[1,1]]))^2&/@MapThread[{#1,#2}&,{data,errors}]]-1)^2,
od>=0&&od<=10*Mean[errors]
},
{od},
Method->{"NelderMead","PostProcess"->False},
AccuracyGoal->2,
PrecisionGoal->2(*,
EvaluationMonitor\[RuleDelayed]Print[od]*)
];

(*SmoothFit=NMinimize[
{
(SquaredErrorFun[p]-Length[data])^2,
p\[GreaterEqual]0&&p\[LessEqual]1
},
{p,0,1(*1/Mean[(1/#^2)&/@errors]*)},
AccuracyGoal\[Rule]1,
PrecisionGoal\[Rule]N[-Log10[Sqrt[2]/Sqrt[Length[data]]]],
EvaluationMonitor\[RuleDelayed]Print[FullForm[p]],
Method\[Rule]{"NelderMead","PostProcess"\[Rule]False}
][[2]];*)

(*SmoothFit=FindRoot[
SquaredErrorFun[p]\[Equal]Length[data],
{p,0,0,1},
AccuracyGoal\[Rule]1,
PrecisionGoal\[Rule]-Log10[Sqrt[2]/Sqrt[Length[data]]](*,
EvaluationMonitor\[RuleDelayed]Print[FullForm[p]]*)
];*)
(*SmoothFit={p\[Rule]1};*)

(*SquaredErrorFun[p]/.SmoothFit//Print;*)

(*SmoothFit=NMinimize[
(p/.SmoothFit)*SquaredErrors+(1-(p/.SmoothFit))*RoughnessFun[deg,knots],
ParamList
];

SmoothFit=BasisFun/.SmoothFit[[2]];*)
{
Evaluate[SmoothFit/.XVar->(#-rescaling[[1]])/(rescaling[[2]]-rescaling[[1]])]&,
od/.ODFit[[2]]
}
)
]


(* ::Input::Initialization:: *)
SSBCalculateCalibrationFit::usage="SSBCalculateCalibrationFit[datatable,n] calculates the calibration fit of the standards in data which was obtained using the standard-sample-bracketing technique. The integer n specifies the number of iterations for the MonteCarloSimulation used to estimate uncertainties.";
SSBCalculateCalibrationFit[data_DataTable,n_]:=Module[
{STDs,Splines,ScaleFactors},
(
STDs=data[[1;;-1;;2]][All,<|#,"Time"->AbsoluteTime[#["Time"]]|>&];
ScaleFactors={Min[STDs[[All,1]]],Max[STDs[[All,1]]]};
QueryTo[STDs,All,<|#,"Time"->(#["Time"]-ScaleFactors[[1]])/(ScaleFactors[[2]]-ScaleFactors[[1]])|>&];

$DistributedContexts=Automatic;
<|
#->Module[{ColDat=NV@STDs[All,{"Time",#}]},
(
Splines=MonteCarloSimulation[
Hold[SmoothSplines[d,errs,3,knots,sf]]/.{
errs:>With[{eval=(DVSEs@ColDat)[[All,2]]},eval/;True],
knots:>With[{eval=ColDat[[All,1]]},eval/;True],
sf:>With[{eval=ScaleFactors},eval/;True]
},
n,
{d->ColDat},
UseParallelKernels->True
];
(*Splines=SmoothSplines[DVMeans@ColDat,(DVSEs@ColDat)[[All,2]],3,ColDat[[All,1]],ScaleFactors];*)

Splines
)
]&/@DeleteCases[ColNames[STDs],"Time"]
|>
)
]


(* ::Input::Initialization:: *)
SSBCalculateDeltas::usage="SSBCalculateDeltas[datatable,standardFit] calculates the \[Delta] values for samples in data which was obtained using the standard-sample-bracketing technique, where the standardFit was previously calculated using SSBCalculateCalibrationFit.";
SSBCalculateDeltas[data_DataTable,standardFit_]:=Module[
{Samples},
(
Samples=data[[2;;-2;;2]][All,<|#,"Time"->AbsoluteTime[#["Time"]]|>&];

Samples[All,Module[{RowV=#},(
<|(
#->((NDV[RowV[[#,1]],RowV[[#,2]]+standardFit[[#]]["SingleMean"][[2]]])/NDV[standardFit[[#]]["SingleMean"][[1]][RowV["Time"]],StandardDeviation@With[{x=RowV["Time"]},#[[1]][x]&/@(standardFit[[#]]["SimulationData"][[All,1]])]]-1)*10^3
)&/@DeleteCases[ColNames[Samples],"Time"]
|>
)]&][All,KeyMap["\[Delta]"<>#<>" / permil"&,#]&]

(*Samples[All,((NDV[#[[2,1]],#[[2,2]]+standardFit["SingleMean"][[2]]])/NDV[standardFit["SingleMean"][[1]][#[[1]]],StandardDeviation@With[{x=#[[1]]},#[[1]][x]&/@(standardFit["SimulationData"][[All,1]])]]-1)*10^3&]*)

(*DTPlot[STDs[[All,{1,3}]],PlotRange\[Rule]All]*)
)
]


(* ::Input::Initialization:: *)
SSBCalculateMeanRatios::usage="SSBCalculateMeanRatios[datatable,{\!\(\*SubscriptBox[\(column\), \(1\)]\)/\!\(\*SubscriptBox[\(column\), \(2\)]\),,\!\(\*SubscriptBox[\(column\), \(3\)]\)/\!\(\*SubscriptBox[\(column\), \(4\)]\),...}] calculates the mean isotope ratios for the standards and samples from the specified columns in data which was obtained using the standard-sample-bracketing technique, by first substracting blank measurements from the standards and samples and filtering outliers.";
SSBCalculateMeanRatios[data_DataTable,ratios_List]:=Module[
{Isotopes,ListRatios,ReplRules,OutDat},
(
Isotopes=DeleteDuplicates@Cases[(List@@#&/@ratios)//Flatten,_String,2];
ListRatios=(Cases[List@@(#),_String,2]//Reverse)&/@ratios;
ReplRules=Association[((#[[1]]<>"/"<>#[[2]])->#[[1]]/#[[2]]&/@ListRatios)];

OutDat=data[All,All,Isotopes];
OutDat=SSBSubstractBlanks[OutDat];
QueryTo[OutDat,All,All,ReplRules/.#&];
QueryTo[OutDat,All,FilterOutliers];
QueryTo[OutDat,All,MeanAndStandardError];

JoinT[data[[2;;-2;;2]][All,All,AbsoluteTime[#["Time"]]&][All,Mean][All,<|"Time"->DateObject[#]|>&],OutDat]
)
]


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
Protect[Evaluate[$Context<>"*"]];
EndPackage[];
