(* ::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`",{
"SciDraw`",
"StyleOptions`",
"CustomTicks`",
"NETLink`",
"ErrorBarPlots`"
}];


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


(* ::Input::Initialization:: *)
IV;
DV;
NDV;
SetConfidenceLevel;

DVEvaluate;
DVE;

DVMeans;
DVSEs;
DVCIs;
NumericDVQ;

MonteCarloSimulation;
MonteCarloResult;
MCEvaluate;
MCE;

EffectiveVarianceFit;

DataTable;
ToArray;

ReplaceDatum;

QueryTo;
KeyFind;
KeyExcept;
GroupByKey;
JoinT;

DefineItemErrors;

FV;
NV;
ColNames;
NK;

FoldValues;
PartitionAssociation;

LSDataStyles;
Fig;
SciDrawPlot;
FigColumn;
FigGrid;
AnchorName;

DTPlot;

VarKey;

MeanAndStandardError;
MeanAndStandardDeviation;
FixedEffectMetaAnalysis;
RandomEffectsMetaAnalysis;

(*MeanAndStandardErrorWithEV;
MeanAndStandardDeviationWithEV;*)

(*EVMultiplyE;*)


(* ::Input::Initialization:: *)
$UsePlusMinusForm;


(* ::Input::Initialization:: *)
CorrelatedErrors;
UseParallelKernels;
AssumeND;
MatchPattern;

XLogScale;
YLogScale;

XErrorBars;
YErrorBars;

FigPanelExtras;
ShowLegend;


(* ::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:: *)
SetAttributes[PrintTiming,HoldFirst];
PrintTiming[f_]:=Module[{res},res=Timing[f];Print[res[[1]]];res[[2]]];


(* ::Input::Initialization:: *)
NumPattern=_Real|_Rational|_Integer;
NumPatternQ[x_]:=MatchQ[x,NumPattern]; 
PosNumPatternQ[x_]:=MatchQ[x,NumPattern]&&x>=0;


(* ::Input::Initialization:: *)
$UsePlusMinusForm::usage="$UsePlusMinusForm is a boolean value specifying whether uncertanties are displayed with the \[PlusMinus] form or with parentheses.";
$UsePlusMinusForm=False;
$ConfidenceLevel=0.95;
$ZValue=InverseCDF[NormalDistribution[],1-(1-$ConfidenceLevel)/2];


(* ::Input::Initialization:: *)
SetConfidenceLevel::usage="SetConfidenceLevel[\!\(\*
StyleBox[\"n\",\nFontSlant->\"Italic\"]\)] sets the confidence level for the DataAnalysis functions, for example to control the output display. The default confidence level is 0.95.";
SetConfidenceLevel[v_:NumPattern]:=(
$ConfidenceLevel=Abs[v];
$ZValue=InverseCDF[NormalDistribution[],1-(1-Abs[v])/2];
v
);


(* ::Input::Initialization:: *)
DV::usage="DV[v,se] represents a value v associated with a standard error se.
DV[v,se,i] represents a value associated with an uncertainty in the form of mean v corresponding and standard error \!\(\*
StyleBox[\"se\",\nFontSlant->\"Italic\"]\) as well as symmetric relative confidence interval i.
DV[v,se,{il,iu}] specifies a confidence interval from v-il to v+iu.";
DV[v_?NumPatternQ,se_?NumPatternQ]:=DV[v,se,Null,Unique[Global`dv]];
DV[v_?NumPatternQ,se_?NumPatternQ,i_?PosNumPatternQ]:=DV[v,se,i,Unique[Global`dv]];
DV[v_?NumPatternQ,se_?NumPatternQ,i:{rl_?PosNumPatternQ,ru_?PosNumPatternQ}]:=DV[v,se,i,Unique[Global`dv]];

DVMeans::usage="DVMeans[expr] replaces all DV, NDV and IV objects in an expression with their mean value.";
DVMeans[expr_]:=ReplaceAllDT[expr,{NDV[v_,se_,n_]:>v,DV[v_,se_,i_,n_]:>v,IV[v_,i_]:>v}];

DVSEs::usage="DVSEs[expr] replaces all DV and NDV objects in an expression with their standard error.";
DVSEs[expr_]:=ReplaceAllDT[expr,{NDV[v_,se_,n_]:>se,DV[v_,se_,i_,n_]:>se}];

DVCIs::usage="DVCIs[expr] gives the confidence intervals for DV and NDV values."
DVCIs[expr_]:=ReplaceAllDT[expr,{
NDV[v_,se_,n_]:>{v-se*$ZValue,v+se*$ZValue},
DV[v_,se_,{rl_?PosNumPatternQ,ru_?PosNumPatternQ},n_]:>{v-rl,v+ru},
DV[v_,se_,i_?PosNumPatternQ,n_]:>{v-i,v+i},
DV[v_,se_,n_]:>{v-se*$ZValue,v+se*$ZValue}
}];

NumericDVQ::usage="NumericDVQ[expr] gives True if expr is a real number or an DV or NDV object with real mean and standard error, and False otherwise.";
NumericDVQ[expr_]:=MatchQ[expr,NumPattern]||
((Head[expr]===NDV||Head[expr]===DV)&&MatchQ[expr[[1]],NumPattern]&&MatchQ[expr[[2]],NumPattern]&&expr[[2]]>0);


(* ::Input::Initialization:: *)
MakeBoxes[dv:DV[v_?NumPatternQ,se_?PosNumPatternQ,Null,_],StandardForm]^:=Interpretation[
Style[SForm[v,se],ShowStringCharacters->False],dv]//ToBoxes;

MakeBoxes[DV[v_?NumPatternQ,se_?PosNumPatternQ,Null,_],OutputForm]^:=OutputSForm[v,se]//ToBoxes;

ToString[DV[v_?NumPatternQ,se_?PosNumPatternQ,Null,_]]^:=OutputSForm[v,se];

MakeBoxes[dv:DV[v_?NumPatternQ,se_?PosNumPatternQ,i_?PosNumPatternQ,_],StandardForm]^:=Interpretation[
Style[If[TrueQ[$UsePlusMinusForm],PMForm[v,i],ParenthesisForm[v,i]],ShowStringCharacters->False],dv]//ToBoxes;

MakeBoxes[DV[v_?NumPatternQ,se_?PosNumPatternQ,i_?PosNumPatternQ,_],OutputForm]^:=If[TrueQ[$UsePlusMinusForm],OutputPMForm[v,i],OutputParenthesisForm[v,i]]//ToBoxes;

ToString[DV[v_?NumPatternQ,se_?PosNumPatternQ,i_?PosNumPatternQ,_]]^:=If[TrueQ[$UsePlusMinusForm],OutputPMForm[v,i],OutputParenthesisForm[v,i]];

MakeBoxes[dv:DV[v_?NumPatternQ,se_?PosNumPatternQ,{rl_?PosNumPatternQ,ru_?PosNumPatternQ},_],StandardForm]^:=Interpretation[
Style[If[TrueQ[$UsePlusMinusForm],PMForm[v,rl,ru],ParenthesisForm[v,rl,ru]],ShowStringCharacters->False],dv]//ToBoxes;

MakeBoxes[DV[v_?NumPatternQ,se_?PosNumPatternQ,{rl_?PosNumPatternQ,ru_?PosNumPatternQ},_],OutputForm]^:=If[TrueQ[$UsePlusMinusForm],OutputPMForm[v,rl,ru],OutputParenthesisForm[v,rl,ru]]//ToBoxes;

ToString[DV[v_?NumPatternQ,se_?PosNumPatternQ,{rl_?PosNumPatternQ,ru_?PosNumPatternQ},_]]^:=If[TrueQ[$UsePlusMinusForm],OutputPMForm[v,rl,ru],OutputParenthesisForm[v,rl,ru]];


(* ::Input::Initialization:: *)
NDV::usage="NDV[v,se] represents a value associated with an uncertainty in the form of a normal distribution with mean v and standard error se.";
NDV[v_?NumPatternQ,se_?PosNumPatternQ]:=NDV[v,se,Unique[Global`dv]];
NDV[v_?NumPatternQ,0]:=v;
NDV[v_?NumPatternQ,0.]:=v;


(* ::Input::Initialization:: *)
MakeBoxes[ndv:NDV[v_?NumPatternQ,se_?PosNumPatternQ,_],StandardForm]^:=Interpretation[
Style[If[TrueQ[$UsePlusMinusForm],PMForm[v,se*$ZValue],ParenthesisForm[v,se*$ZValue]],ShowStringCharacters->False],ndv]//ToBoxes;

MakeBoxes[NDV[v_?NumPatternQ,se_?PosNumPatternQ,_],OutputForm]^:=If[TrueQ[$UsePlusMinusForm],OutputPMForm[v,se*$ZValue],OutputParenthesisForm[v,se*$ZValue]]//ToBoxes;

ToString[NDV[v_?NumPatternQ,se_?PosNumPatternQ,_]]^:=If[TrueQ[$UsePlusMinusForm],OutputPMForm[v,se*$ZValue],OutputParenthesisForm[v,se*$ZValue]];


(* ::Input::Initialization:: *)
IV::usage="IV[v,r] represents a value associated with a confidence interval that ranges from v-r to v+r.
IV[v,{r1,r2}] represents a value associated with a confidence interval from v-r1 to v+r2.
IV[NDV[..]] gives the IV for a NDV object.";
IV[NDV[v_?NumPatternQ,se_?PosNumPatternQ,_]]:=IV[v,se*$ZValue];


(* ::Input::Initialization:: *)
OutputParenthesisForm[v_,r_]:=Module[
{PF},
(
PF=ParenthesisForm[v,r];
PF=PF/.{Superscript[10,Pow_]:>("10^"<>ToString[Pow])};
ToString[PF]
)
];

ParenthesisForm[v_,r_]:=Module[{Pow,n=Ceiling[-Log10[r]]+1},
(
If[v!=0&&(Abs[v]<10^-3||Abs[v]>=10^4),
Pow=Ceiling[-Log10[Abs[v]]];
n=Ceiling[-Log10[r*10^Pow]]+1;

If[n<=0,
Row[{Round[v*10^Pow],"(",Round[r*10^Pow],")\[Times]",Superscript[10,-Pow]}],
If[!NumberQ[n],
v,
Row[{NumberForm[v*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"(",Round[r*10^Pow*10^n],")\[Times]",Superscript[10,-Pow]}]
]
],
If[TrueQ[n<=0],
Row[{Round[v],"(",Round[r],")"}],
If[!NumberQ[n],
v,
Row[{NumberForm[v//N,{\[Infinity],n},ExponentFunction->(Null&)],"(",Round[r*10^n],")"}]
]
]
]
)
];

OutputParenthesisForm[v_,rl_,ru_]:=Module[
{PF},
(
PF=ParenthesisForm[v,rl,ru];
PF=PF/.{Superscript[10,Pow_]:>("10^"<>ToString[Pow])};
ToString[PF]
)
];

ParenthesisForm[v_,rl_,ru_]:=Module[{Pow,n=Ceiling[-Log10[Min[rl,ru]]]+1},
(
If[v!=0&&(Abs[v]<10^-3||Abs[v]>=10^4),
Pow=Ceiling[-Log10[Abs[v]]];
n=Ceiling[-Log10[Min[rl,ru]*10^Pow]]+1;

If[TrueQ[n<=0],
Row[{Round[v*10^Pow],"(-",Round[rl*10^Pow],"+",Round[ru*10^Pow],")\[Times]",Superscript[10,-Pow]}],
If[!NumberQ[n],
v,
Row[{NumberForm[v*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"(-",Round[rl*10^Pow*10^n],"+",Round[ru*10^Pow*10^n],")\[Times]",Superscript[10,-Pow]}]
]
],
If[TrueQ[n<=0],
Row[{Round[v],"(-",Round[rl],"+",Round[ru],")"}],
If[!NumberQ[n],
v,
Row[{NumberForm[v//N,{\[Infinity],n},ExponentFunction->(Null&)],"(-",Round[rl*10^n],"+",Round[ru*10^n],")"}]
]
]
]
)
]


(* ::Input::Initialization:: *)
OutputPMForm[v_,r_]:=Module[
{PF},
(
PF=PMForm[v,r];
PF=PF/.{Superscript[10,Pow_]:>("10^"<>ToString[Pow])};
ToString[PF]
)
];

PMForm[v_,r_]:=Module[{Pow,n=Ceiling[-Log10[r]]+1},
(
If[v!=0&&(Abs[v]<10^-3||Abs[v]>=10^4),
Pow=Ceiling[-Log10[Abs[v]]];
n=Ceiling[-Log10[r*10^Pow]]+1;

If[n<=0,
Row[{Round[v*10^Pow],"\[PlusMinus]",Round[r*10^Pow],"\[Times]",Superscript[10,-Pow]}],
If[!NumberQ[n],
v,
Row[{NumberForm[v*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"\[PlusMinus]",NumberForm[r*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"\[Times]",Superscript[10,-Pow]}]
]
],
If[TrueQ[n<=0],
Row[{Round[v],"\[PlusMinus]",Round[r]}],
If[!NumberQ[n],
v,
Row[{NumberForm[v//N,{\[Infinity],n},ExponentFunction->(Null&)],"\[PlusMinus]",NumberForm[r//N,{\[Infinity],n},ExponentFunction->(Null&)]}]
]
]
]
)
];

OutputPMForm[v_,rl_,ru_]:=Module[
{PF},
(
PF=PMForm[v,rl,ru];
PF=PF/.{Superscript[10,Pow_]:>("10^"<>ToString[Pow])};
ToString[PF]
)
];

PMForm[v_,rl_,ru_]:=Module[{Pow,n=Ceiling[-Log10[Min[rl,ru]]]+1},
(
If[v!=0&&(Abs[v]<10^-3||Abs[v]>=10^4),
Pow=Ceiling[-Log10[Abs[v]]];
n=Ceiling[-Log10[Min[rl,ru]*10^Pow]]+1;

If[n<=0,
Row[{Round[v*10^Pow],"\[PlusMinus][-",Round[rl*10^Pow],"+",Round[ru*10^Pow],"]\[Times]",Superscript[10,-Pow]}],
If[!NumberQ[n],
v,
Row[{NumberForm[v*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"\[PlusMinus][-",NumberForm[rl*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"+",NumberForm[ru*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"]\[Times]",Superscript[10,-Pow]}]
]
],
If[TrueQ[n<=0],
Row[{Round[v],"\[PlusMinus][-",Round[rl],"+",Round[ru],"]"}],
If[!NumberQ[n],
v,
Row[{NumberForm[v//N,{\[Infinity],n},ExponentFunction->(Null&)],"\[PlusMinus][-",NumberForm[rl//N,{\[Infinity],n},ExponentFunction->(Null&)],"+",NumberForm[ru//N,{\[Infinity],n},ExponentFunction->(Null&)],"]"}]
]
]
]
)
]


(* ::Input::Initialization:: *)
OutputSForm[v_,se_]:=Module[
{PF},
(
PF=PMForm[v,se];
PF=PF/.{Superscript[10,Pow_]:>("10^"<>ToString[Pow])};
ToString[PF]
)
];

SForm[v_,se_]:=Module[{Pow,n=Ceiling[-Log10[se]]+1},
(
If[v!=0&&(Abs[v]<10^-3||Abs[v]>=10^4),
Pow=Ceiling[-Log10[Abs[v]]];
n=Ceiling[-Log10[se*10^Pow]]+1;

If[n<=0,
Row[{Round[v*10^Pow],"(SE=",Round[se*10^Pow],"\[Times]",Superscript[10,-Pow],")"}],
If[!NumberQ[n],
v,
Row[{NumberForm[v*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"(SE=",NumberForm[se*10^Pow//N,{\[Infinity],n},ExponentFunction->(Null&)],"\[Times]",Superscript[10,-Pow],")"}]
]
],
If[TrueQ[n<=0],
Row[{Round[v],"(SE=",Round[se],")"}],
If[!NumberQ[n],
v,
Row[{NumberForm[v//N,{\[Infinity],n},ExponentFunction->(Null&)],"(SE=",NumberForm[se//N,{\[Infinity],n},ExponentFunction->(Null&)],")"}]
]
]
]
)
]


(* ::Input::Initialization:: *)
Options[DVEvaluate]={CorrelatedErrors->False,AssumeND->True};
CorrelatedErrors::usage="CorrelatedErrors is an option for DVEvaluate to specify whether to assume that all random errors are perfectly correlated (True) or completely uncorrelated (False).";
AssumeND::usage="AssumeND is an option for DVEvaluate to specify whether to assume that the resulting distribution is a normal distribution (True) or not (False).";
(*DVEvaluate[expr_,opts:OptionsPattern[]]:=expr/.{p_?EVExpressionQ\[RuleDelayed]DVEvaluateSingle[p,opts]};*)
(*Attributes[DVEvaluate]={HoldFirst};*)

DVEvaluate::usage="DVEvaluate[expr] evaluates expr while propagating standard errors of DV and NDV values by linearizing the functions in the expression.";
DVEvaluate[expr_,opts:OptionsPattern[]]:=Module[
{Params,Rules,V,S,F},
(
F=expr(*/.EV[v_,s_]\[RuleDelayed]EV[v,s,Unique[Global`ev]]*);
(*F=Hold[expr]/.ev_EV\[RuleDelayed]With[{eval=Evaluate[ev]},eval/;True];*)
Params=Cases[F,_DV|_NDV,{0,Infinity}];

If[Length[Params]>0,
(
Params=DeleteDuplicates[Params];
V=F(*//ReleaseHold*)//DVMeans//Quiet;

If[(*(Quiet[Depth[F(*//ReleaseHold*)]]>2)*)(Quiet[!MatchQ[F,_DV|_NDV]])&&NumericQ[V],
(
Rules=Rule[Unique[dvslot],#]&/@Params;
F=Quiet[F/.(Reverse/@Rules)(*//ReleaseHold*)];

S=If[TrueQ[OptionValue[CorrelatedErrors]],
(
Sqrt[
Plus@@Map[(D[F,#[[1]]]/.DVMeans[Rules])^2*DVSEs[#[[2]]]^2&,Rules]+
2*Plus@@(((D[F,#[[1,1]]]/.DVMeans[Rules])*
(D[F,#[[2,1]]]/.DVMeans[Rules])*
DVSEs[((#[[1,2]])*(#[[2,2]]))]

)&/@Subsets[Cases[Rules,HoldPattern[_->(DV[___]|NDV[___])]],{2}])
]
),
(
Sqrt[Plus@@Map[(D[F,#[[1]]]/.DVMeans[Rules])^2*DVSEs[#[[2]]]^2&,Rules]]
)
]//Quiet;

If[!NumericQ[S],
Map[DVEvaluate[#,opts]&,F](*Replace[Quiet[expr],x_\[RuleDelayed]DVEvaluate[x,opts],{1}]*),
If[S==0,V,If[TrueQ[OptionValue[AssumeND]],NDV[V//N,S//N],DV[V//N,S//N]]]
]//Quiet
),
(
(*Map["fubar"&,F]//Print;Abort[];*)
Map[DVEvaluate[#,opts]&,F]
(*Replace[Quiet[F(*//ReleaseHold*)],x_\[RuleDelayed]DVEvaluate[x,opts],{1}]*)
)
]
)
,
expr
]
)
]

Options[DVE]={CorrelatedErrors->False};
(*Attributes[DVE]={HoldFirst};*)

DVE::usage="DVE[expr] represents a short form of DVEvaluate[expr].";
DVE[expr_,opts:OptionsPattern[]]:=(DVEvaluate[expr,opts]);

(*Options[DVEvaluateSingle]={CorrelatedErrors\[Rule]False};
DVEvaluateSingle[expr_,opts:OptionsPattern[]]:=Module[
{Params,Rules},
(
Params=Cases[expr,_EV,{0,Infinity}];

If[Length[Params]>0,
(
Params=Params/.EV[v_,s_]\[RuleDelayed]EV[v,s,Unique[Global`ev]];
Params=DeleteDuplicates[Params];

Rules=Rule[#,Unique[dvslot]]&/@Params;
expr/.Rules/.If[TrueQ[OptionValue[CorrelatedErrors]],CorrelatedEVRules[Sequence@@(Reverse/@Rules)],EVRules[Sequence@@(Reverse/@Rules)]]
)
,
expr
]
)
]*)


(* ::Input::Initialization:: *)
$DistributionList=ToExpression/@Names["*Distribution"];
$DistributionRandomRules=((ToExpression["x_"<>ToString[#]]):>RandomVariate[x]&/@$DistributionList);
$DistributionMeanRules=((ToExpression["x_"<>ToString[#]]):>Mean[x]&/@$DistributionList);


(* ::Input::Initialization:: *)
RandomData[dat_]:=Module[
{},
(
ReplaceAllDT[ReplaceAllDT[dat,$DistributionRandomRules],(NDV[v_,se_,_]:>If[se>0,RandomVariate[NormalDistribution[v,se]],v])]
)
];


(* ::Input::Initialization:: *)
ReplaceMeans[dat_]:=Module[
{},
(
DVMeans[ReplaceAllDT[dat,$DistributionMeanRules]]
)
];


(* ::Input::Initialization:: *)
ReleaseHoldAt[expr_,partspec_]:=ReplacePart[expr,partspec->Extract[expr,partspec]];
EvaluateAssocValues[expr_]:=ReplaceAll[expr,a:Association[___]:>((EvaluateAssocValues[#]&/@a))];
ReplaceAllDT[expr_,rules_]:=(expr/.rules)/.d_DataTable:>(d/.rules);


(* ::Input::Initialization:: *)
(*SetAttributes[MonteCarloSimulation,HoldFirst];*)
Options[MonteCarloSimulation]={UseParallelKernels->False,MatchPattern->Blank[],MissingBehavior->Automatic};
UseParallelKernels::usage="UseParallelKernels is an option for MonteCarloSimulation[...] which specifies whether parallel kernels should be used for the simulation.";
MatchPattern::usage="MatchPattern is an option for MonteCarloSimulation[...] which specifies that only function results matching this pattern should be used (default is Blank[]).";

MonteCarloSimulation::minnum="The number of simulations i has to be greater than 1.";

MonteCarloSimulation::usage="MonteCarloSimulation[expr,i,{\!\(\*SubscriptBox[\(param\), \(1\)]\)\[Rule]\!\(\*SubscriptBox[\(distribution\), \(1\)]\),\!\(\*SubscriptBox[\(param\), \(2\)]\)\[Rule]\!\(\*SubscriptBox[\(distribution\), \(2\)]\)}] carries out a Monte Carlo simulation with i iterations, to propagate uncertainties (in NDV objects or standard Mathematica Distribution objects) in expr. Parameters in expr may be separately specified as rules \!\(\*SubscriptBox[\(param\), \(i\)]\)\[Rule]\!\(\*SubscriptBox[\(distribution\), \(i\)]\), where \!\(\*SubscriptBox[\(distribution\), \(i\)]\) can be either an NDV object or a standard Mathematica distribution object.";
MonteCarloSimulation[f_,i_Integer,params_,opts:OptionsPattern[]]:=Module[
{TmpRes,Counter,StartCounter,SingleRes,MP=OptionValue[MatchPattern]},
(
If[MissingBehavior===Automatic&&(!FreeQ[params,Missing[___]]||!FreeQ[Hold[f],Missing[___]]),Return[Missing[]]];

If[i<2,ErrorMessage[MonteCarloSimulation::minnum];];

If[TrueQ[OptionValue[UseParallelKernels]],
(
Counter=0;
StartCounter=0;
SetSharedVariable[Counter,StartCounter];

Monitor[TmpRes=Flatten[ParallelEvaluate[Module[{Res={},ThisRes,P},(
While[StartCounter<i,
StartCounter++;
P=EvaluateAssocValues@(params//RandomData);

ThisRes=(*EvaluateAssocValues@*)((f//RandomData)/.P//ReleaseHold);
If[MatchQ[Evaluate[ThisRes],MP],
(
AppendTo[Res,{ThisRes,P}];
Counter++;
),(
StartCounter--;
)];
];
Res
)],DistributedContexts:>If[ListQ[$DistributedContexts],Append[$DistributedContexts,"DataAnalysis`Private`"],If[StringQ[$DistributedContexts],{$DistributedContexts,"DataAnalysis`Private`"},$DistributedContexts]]],1];,Column[{
ToString[Counter]<>"/"<>ToString[i],
ProgressIndicator[Counter,{0,i}]
}]];

If[Length[TmpRes]>i,TmpRes=TmpRes[[;;i]]];
),
(
Counter=0;
TmpRes={};
Monitor[Module[{ThisRes,P},(
While[Counter<i,
P=EvaluateAssocValues@(params//RandomData);
ThisRes=(*EvaluateAssocValues@*)((f//RandomData)/.P//ReleaseHold);
If[MatchQ[ThisRes,MP],
(
AppendTo[TmpRes,{ThisRes,P}];
Counter++;
)];
];
)],
Column[{
ToString[Counter]<>"/"<>ToString[i],
ProgressIndicator[Counter,{0,i}]
}]
];
)
];

SingleRes=(f//ReplaceMeans)/.(params//ReplaceMeans//EvaluateAssocValues)//ReleaseHold(*//EvaluateAssocValues*);

MonteCarloResult[
{SingleRes,TmpRes}
(*If[AssociationQ[SingleRes],{Map[#&,SingleRes,Infinity],Map[#&,TmpRes,{3}]},{SingleRes,TmpRes}]*)
]
)
]


(* ::Input::Initialization:: *)
MakeBoxes[res:MonteCarloResult[data_List],StandardForm]^:=MonteCarloResult[Panel[Row[{"n = ",data[[2]]//Length}],FrameMargins->5]]//ToBoxes;

MonteCarloResult::usage="MonteCarloResult[...] represents the symbolic result obtained from a MonteCarloSimulation.";
MonteCarloResult::unknownproperty="The property \"`1`\" is not defined.";
MonteCarloResult::resultform="This output can only be generated if the result of the function is numeric, a flat list, or a flat association.";
MonteCarloResult::reslist="This output can only be generated for numerical results (e.g. not for lists or associations).";
MonteCarloResult::paramlist="Sensitifity plots can only be generated for numerical parameters (e.g. not if one parameter is a list or association).";
MonteCarloResult::notnormal="A DistributionFitTest indicates that the simulated results are not normally distributed. This NDV value has to be treated with caution, as NDV value are assumed to be normally distributed.";
MonteCarloResult[data_][property_String]:=Module[
{},
(
Switch[property,
"SingleMean",
data[[1]],

"SingleMeanDV",
CheckResultForm[data[[1]]];
Switch[data[[1]],
_List,MapThread[(DV[#1,StandardDeviation[#2],{-1,1}*(QuantileCI[#2]-#1)])&,{data[[1]],data[[2,All,1]]//Transpose},Depth[data[[1]]]-1],
_Association,MapAssoc[(DV[#1,StandardDeviation[#2],{-1,1}*(QuantileCI[#2]-#1)])&,data],
NumPattern,(DV[data[[1]],StandardDeviation[data[[2,All,1]]],{-1,1}*(QuantileCI[data[[2,All,1]]]-data[[1]])])
],

"SingleMeanNDV",
CheckResultForm[data[[1]]];
Switch[data[[1]],
_List,MapThread[(CheckNormality[#2];NDV[#1,StandardDeviation[#2]])&,{data[[1]],data[[2,All,1]]//Transpose},Depth[data[[1]]]-1],
_Association,MapAssoc[(CheckNormality[#2];NDV[#1,StandardDeviation[#2]])&,data],
NumPattern,(CheckNormality[data[[2,All,1]]];NDV[data[[1]],StandardDeviation[data[[2,All,1]]]])
],

"QuantileCI",
CheckResultForm[data[[1]]];
Switch[data[[1]],
_List,MapThread[(QuantileCI[#2])&,{data[[1]],data[[2,All,1]]//Transpose},Depth[data[[1]]]-1],
_Association,MapAssoc[(QuantileCI[#2])&,data],
NumPattern,(QuantileCI[data[[2,All,1]]])
],

"DistributionMean",
Mean[data[[2,All,1]]],

"DistributionMeanDV",
CheckResultForm[data[[2,1,1]]];
Switch[data[[2,1,1]],
_List,DV[Mean[#],StandardDeviation[#],{-1,1}*(QuantileCI[#]-Mean[#])]&/@(data[[2,All,1]]//Transpose),
_Association,MapAssoc[DV[Mean[#2],StandardDeviation[#2],{-1,1}*(QuantileCI[#2]-Mean[#2])]&,data],
NumPattern,DV[Mean[data[[2,All,1]]],StandardDeviation[data[[2,All,1]]],{-1,1}*(QuantileCI[data[[2,All,1]]]-Mean[data[[2,All,1]]])]
],

"DistributionMeanNDV",
CheckResultForm[data[[2,1,1]]];
Switch[data[[2,1,1]],
_List,(CheckNormality[#];NDV[Mean[#],StandardDeviation[#]])&/@(data[[2,All,1]]//Transpose),
_Association,MapAssoc[(CheckNormality[#2];NDV[Mean[#2],StandardDeviation[#2]])&,data],
NumPattern,(CheckNormality[data[[2,All,1]]];NDV[Mean[data[[2,All,1]]],StandardDeviation[data[[2,All,1]]]])
],

"FindDistribution",
CheckResultForm[data[[2,1,1]]];
Switch[data[[2,1,1]],
_List,FindDistribution[#]&/@(data[[2,All,1]]//Transpose),
_Association,MapAssoc[(FindDistribution[#2])&,data],
NumPattern,FindDistribution[data[[2,All,1]]]
],

"EmpiricalDistribution",
CheckResultForm[data[[2,1,1]]];
Switch[data[[2,1,1]],
_List,EmpiricalDistribution[#]&/@(data[[2,All,1]]//Transpose),
_Association,MapAssoc[EmpiricalDistribution[#2]&,data],
NumPattern,EmpiricalDistribution[data[[2,All,1]]]
],

"StandardDeviation",
CheckResultForm[data[[2,1,1]]];
Switch[data[[2,1,1]],
_List,StandardDeviation[data[[2,All,1]]],
_Association,MapAssoc[StandardDeviation[#2]&,data],
NumPattern,StandardDeviation[data[[2,All,1]]]
],

"Histogram",
CheckResultForm[data[[2,1,1]]];
Switch[data[[2,1,1]],
_List,Histogram[#]&/@(data[[2,All,1]]//Transpose),
_Association,MapAssoc[Histogram[#2]&,data],
NumPattern,Histogram[data[[2,All,1]]]
],

"SensitivityPlots",Module[
{Params},
(
If[!NumericQ[data[[1]]],ErrorMessage[MonteCarloResult::reslist]];
If[!(And@@(NumericQ/@data[[2,1,2,All,2]])),ErrorMessage[MonteCarloResult::paramlist]];

Params=data[[2,1,2,All,1]];
Multicolumn[#,1,Spacings->5]&@(ListPlot[
(MapThread[{#1,#2}&,{#/.data[[2,All,2]],data[[2,All,1]]}]),
Frame->True,
FrameStyle->Black,
PlotRange->All,
ImageSize->450,
PlotMarkers->$DefaultPlotMarkers,
PlotStyle->$DefaultPlotColors,
FrameTicksStyle->12,
LabelStyle->12,
Axes->False,
PlotLabel->Style[#,16,Black]
]&/@Params)
)
],

"SimulationData",
data[[2]],

_,
ErrorMessage[MonteCarloResult::unknownproperty,property]
]
)
]

MapAssoc[f_,data_]:=Association@(Module[{KeyName,Res},
(
KeyName=#;
Res=#[KeyName]&/@data[[2,All,1]];

KeyName->If[AssociationQ[data[[1]][KeyName]],
MapSubAssoc[f,data[[1]][KeyName],Res],
f[data[[1]][KeyName],Res]
]
)
]&/@Keys[data[[1]]]);

MapSubAssoc[f_,data0_,data_]:=Association@(Module[{KeyName,Res},
(
KeyName=#;
Res=#[KeyName]&/@data;

KeyName->If[AssociationQ[data0[KeyName]],
MapSubAssoc[f,data0[KeyName],Res],
f[data0[KeyName],Res]
]
)
]&/@Keys[data0]);

CheckResultForm[res_]:=If[MatchQ[res,(_Real|_Integer|_Rational)|(<|(_->_Real|_Integer|_Rational|_Association)...|>)|({(_Real|_Integer|_Rational)...})(*(x_/;(ListQ[x]&&ArrayQ[x]))*)],True,ErrorMessage[MonteCarloResult::resultform]];

CheckNormality[data_]:=If[DistributionFitTest[data]<(1-$ConfidenceLevel),Message[MonteCarloResult::notnormal]];

QuantileCI[data_]:=Quantile[data,{(1-$ConfidenceLevel)/2,1-(1-$ConfidenceLevel)/2}];

QuantileMeanRelCI[data_]:=Abs[Quantile[data,{(1-$ConfidenceLevel)/2,1-(1-$ConfidenceLevel)/2}]-Mean[data]];
(*Module[{BootSample},
(
BootSample=Table[Mean[RandomChoice[data,Length@data]],{n}];
Abs[Quantile[BootSample,{(1-$ConfidenceLevel)/2,1-(1-$ConfidenceLevel)/2}]-Mean[BootSample]]
)]*)


(* ::Input::Initialization:: *)
Options[MCEvaluate]=Options[MonteCarloSimulation];
(*SetAttributes[MCEvaluate,HoldFirst];*)

MCEvaluate::usage="MCEvaluate[expr] evaluates expr while propagating errors in NDV values using a Monte Carlo simulation with 1000 iterations.
MCEvaluate[expr,i] evaluates expr while propagating errors in NDV values using a Monte Carlo simulation with i iterations.";
MCEvaluate[expr_,i_:10000,opts:OptionsPattern[]]:=MonteCarloSimulation[expr,i,{},opts]["SingleMeanDV"];

Options[MCE]=Options[MCEvaluate];
(*SetAttributes[MCE,HoldFirst];*)
MCE[expr_,i_:1000,opts:OptionsPattern[]]:=MCEvaluate[expr,i,opts];


(* ::Input::Initialization:: *)
Options[EffectiveVarianceFit]=Options[NonlinearModelFit];

EffectiveVarianceFit::usage="EffectiveVarianceFit[{{\!\(\*SubscriptBox[\(x\), \(11\)]\),\!\(\*SubscriptBox[\(x\), \(12\)]\),\[Ellipsis],\!\(\*SubscriptBox[\(y\), \(1\)]\)},{\!\(\*SubscriptBox[\(x\), \(21\)]\),\!\(\*SubscriptBox[\(x\), \(22\)]\),\[Ellipsis],\!\(\*SubscriptBox[\(y\), \(2\)]\)},\[Ellipsis]},{{s\!\(\*SubscriptBox[\(x\), \(11\)]\),s\!\(\*SubscriptBox[\(x\), \(12\)]\),\[Ellipsis],s\!\(\*SubscriptBox[\(y\), \(1\)]\)},{s\!\(\*SubscriptBox[\(x\), \(21\)]\),s\!\(\*SubscriptBox[\(x\), \(22\)]\),\[Ellipsis],s\!\(\*SubscriptBox[\(y\), \(2\)]\)},\[Ellipsis]},form,{\!\(\*SubscriptBox[\(\[Beta]\), \(1\)]\),\[Ellipsis]},{\!\(\*SubscriptBox[\(x\), \(1\)]\),\[Ellipsis]},i] constructs a nonlinear model with structure form where form depends on variables \!\(\*SubscriptBox[\(x\), \(i\)]\) using the parameters \!\(\*SubscriptBox[\(\[Beta]\), \(i\)]\). An effective variance technique is used to fit the data using errors \!\(\*SubscriptBox[\(sx\), \(ij\)]\) and \!\(\*SubscriptBox[\(sy\), \(j\)]\) in the data points \!\(\*SubscriptBox[\(x\), \(ij\)]\) and \!\(\*SubscriptBox[\(y\), \(j\)]\) using i iterations.
EffectiveVarianceFit[datatable,form,{\!\(\*SubscriptBox[\(\[Beta]\), \(1\)]\),\[Ellipsis]},{\!\(\*SubscriptBox[\(x\), \(1\)]\),\[Ellipsis]},i] constructs a nonlinear model with structure form where form depends on variables \!\(\*SubscriptBox[\(x\), \(i\)]\) using the parameters \!\(\*SubscriptBox[\(\[Beta]\), \(i\)]\). An effective variance technique is used to fit the values and errors in the DataTable using i iterations.
";
EffectiveVarianceFit[data_List,errors_List,form_,params_,vars_,runs_,opts:OptionsPattern[]]:=Module[
{YErrors,FitRes,DSquares},
(
YErrors=errors[[All,-1]];

Do[(
FitRes=NonlinearModelFit[data,form,params,vars,opts,VarianceEstimatorFunction->(1&),Weights->(1/#^2&/@YErrors)];

DSquares=(D[If[ListQ@form,form[[1]],form]/.FitRes["BestFitParameters"],#]^2)&/@vars;
DSquares=((#*DSquares)&/@(errors[[All,;;-2]]^2));

DSquares=MapThread[#1/.#2&,{
DSquares,
MapThread[#1->#2&,{vars,#//Most}]&/@data
}];

YErrors=MapThread[Sqrt[#1+#2]&,{errors[[All,-1]]^2,(Plus@@#&/@DSquares)}];
),{runs}];

FitRes
)
];

EffectiveVarianceFit::noerror="Some of the input data are not NDV items.";
EffectiveVarianceFit[tbl_DataTable,form_,params_,vars_,runs_,opts:OptionsPattern[]]:=Module[
{Dat,VarMap,FitRes},(
Dat=tbl[All,Append[vars,Last@ColNames@tbl]];
Dat=Dat[Select[#,MatchQ[#,<|(_->_?NumericDVQ)..|>]&]&];

VarMap=#->VarKey[#]&/@vars;
Unprotect[DataAnalysis`VarKey];
FitRes=EffectiveVarianceFit[NV@Dat[All,All,DVMeans],NV@Dat[All,All,If[Head[#]===NDV,DVSEs@#,ErrorMessage[EffectiveVarianceFit::noerror]]&],form/.VarMap,params,vars/.VarMap,runs,opts];
FitRes
)];


(* ::Input::Initialization:: *)
(*(* Addition *)
Plus[EV[v1_,s1_],EV[v2_,s2_]]^:=
EV[
v1+v2,
Sqrt[s1^2+s2^2]//N
];
Plus[v1_?(MatchQ[N[#],_Integer|_Rational|_Real]&),EV[v2_,s2_,n2_]]^:=
EV[
v1+v2,
s2
];

(* Multiplication *)
Times[EV[v1_,s1_],EV[v2_,s2_]]^:=
EV[
v1*v2,
Abs[v1*v2]*Sqrt[(s1/v1)^2+(s2/v2)^2]
];
Times[v1_?(MatchQ[N[#],_Integer|_Rational|_Real]&),EV[v2_,s2_,n2_]]^:=
EV[
v1*v2,
N[Abs[v1]*s2]
];

(* Power *)
Power[EV[v1_,s1_],EV[v2_,s2_]]^:=
EV[
v1^v2,
N[Abs[v1^v2]*Sqrt[(v2/v1*s1)^2+(Log[v1]*s2)^2]]
];
Power[EV[v1_,s1_],c_?(MatchQ[N[#],_Integer|_Rational|_Real]&)]^:=
EV[
v1^c,
N[Abs[v1^c*c*s1/v1]]
];
Power[c_?(MatchQ[N[#],_Integer|_Rational|_Real]&),EV[v1_,s1_]]^:=
EV[
c^v1,
N[Abs[c^v1*Log[c]*s1]]
];

(* Log *)
Log[EV[v1_,s1_]]^:=
EV[
Log[v1],
N[Abs[s1/v1]]
];

(* Basic trigonometric functions *)
Sin[EV[v1_,s1_]]^:=EV[Sin[v1],N[Abs[s1*Cos[v1]]]];
Cos[EV[v1_,s1_]]^:=EV[Cos[v1],N[Abs[s1*Sin[v1]]]];
Tan[EV[v1_,s1_]]^:=EV[Tan[v1],N[s1*Abs[Sec[v1]]^2]];*)


(* ::Input::Initialization:: *)
(* Assignments *)

(* Rule *)
(*Rule[lhs_,EV[v_,s_]]^:=Rule[lhs,EV[v,s,Unique[Global`ev]]];*)

(* Replace *)
(*ReplaceAll[expr_,EVRules[rules__]]^:=
Module[{ValueRuleList,Err},
(
ValueRuleList=(List@rules)/.EV[v_,s_,n_]\[RuleDelayed]v;

Err=Sqrt[Plus@@Map[
(D[expr,#[[1]]]/.ValueRuleList)^2*(#[[2]]/.EV[v_,s_,n_]\[RuleDelayed]s)^2&,
Cases[List@rules,HoldPattern[_\[Rule]EV[___]]]
]]//N;

If[Err\[Equal]0,
expr/.ValueRuleList,
EV[
expr/.ValueRuleList,
Err
]
]
)
];
ReplaceAll[expr_,CorrelatedEVRules[rules__]]^:=
Module[{ValueRuleList,Err},
(
ValueRuleList=(List@rules)/.EV[v_,s_,n_]\[RuleDelayed]v;

Err=Sqrt[Plus@@Map[
(D[expr,#[[1]]]/.ValueRuleList)^2*(#[[2]]/.EV[v_,s_,n_]\[RuleDelayed]s)^2&,
Cases[List@rules,HoldPattern[_\[Rule]EV[___]]]
]+
2*Plus@@(((D[expr,#[[1,1]]]/.ValueRuleList)*
(D[expr,#[[2,1]]]/.ValueRuleList)*
((#[[1,2]]/.ValueRuleList)*(#[[2,2]]/.ValueRuleList))/.EV[v_,s_,n_]\[RuleDelayed]s

)&/@Subsets[Cases[List@rules,HoldPattern[_\[Rule]EV[___]]],{2}])
]//N;

If[Err\[Equal]0,
expr/.ValueRuleList,
EV[
expr/.ValueRuleList,
Err
]
]
)
];*)

(* Function *)
(*Function[f_][pre___,EV[v1_,s1_],suff___]^:=Module[
{Args,ValueList,OutF,EVs,Expr,Rules},
(
Args=Join[List@pre,{EV[v1,s1]},List@suff];
Args=Args/.EV[v_,s_]\[RuleDelayed]EV[v,s,Unique[Global`ev]];
(*ValueArgList=Args/.EV[v_,s_]\[RuleDelayed]v;
IDList=Table[Unique[ev],{Count[Args,_EV]}];*)
(*ValueList=Args/.EV[v_,s_,n_]\[RuleDelayed]v;*)

OutF=Function[f][Sequence@@Args];
EVs=Cases[Args,_EV];
Rules=Rule[#,Unique[dvslot]]&/@Cases[Args,_EV];
OutF/.Rules/.EVRules[Sequence@@(Reverse/@Rules)]

(*EV[
Function[f][Sequence@@ValueList],
Sqrt[Plus@@Map[
(
SlotN=#;
(*(Function[D[f,Slot[#]]][Sequence@@ValueList])^2*(Args[[#]]/.EV[v_,s_,n_]\[RuleDelayed]s)^2*)
D[f,Slot[Slot[1]]]
)&,
Position[Args,_EV]//Flatten
]]
]*)

(*Function[f][Sequence@@Args]*)
)
];*)


(* ::Input::Initialization:: *)
(* Addition *)
(*Plus[EV[v1_,s1_,n1_],EV[v1_,s1_,n1_]]^:=
EV[
2*v1,
2*s1,
n1
];
Plus[v1_?(MatchQ[N[#],_Integer|_Rational|_Real]&),EV[v2_,s2_,n2_]]^:=
EV[
v1+v2,
s2,
n2
];*)

(* Multiplication *)
(*Times[EV[v1_,s1_,n1_],EV[v1_,s1_,n1_]]^:=
EV[
v1^2,
Abs[v1^2*2*s1/v1],
n1
];*)

(* Power *)
(*Power[EV[v1_,s1_,n1_],EV[v1_,s1_,n1_]]^:=
EV[
v1^v2,
Abs[v1^v2]*Sqrt[(v2/v1*s1)^2+(Log[v1]*s2)^2],
n1
];*)


(* ::Input::Initialization:: *)
(*EVExpressionQ[expr_]:=(Depth[expr]>2)&&(And@@(MatchQ[#,_Integer|_Real|EV[__]|_?((StringMatchQ[ToString[#],"ev$"~~NumberString])&)|_Plus|_Times]&/@(Level[expr,Infinity])))*)

(*EVExpressionQ[expr_]:=(Depth[expr]>1)&&(
NumericQ[expr/.EV[v_,s_,n_]\[RuleDelayed]v]
)*)


(* ::Input::Initialization:: *)
DataTable::usage="DataTable[data] represents a structured data table, based on and extending the Dataset structure.";
DataTable[data_Association]:=DataTable[Dataset[data]];

DataListQ[tbl_]:=ArrayQ[tbl,2]&&MatchQ[tbl,{{__}..}]

DataTable::wrongdimensions="The list must a full 2-dimensional array.";
DataTable[data_List]:=Module[{Columns,ID,IDs,Out,FFun,CheckQ,i},
(
If[!DataListQ[data],ErrorMessage[DataTable::wrongdimensions]];

Out=data/.(Global`\[Placeholder]->Missing[]);
Columns=Out//First//Rest;

IDs={};
Out=(
ID=#//First;
i=2;
While[MemberQ[IDs,ID],ID=ToString[ID]<>"-"<>ToString[i++]];
AppendTo[IDs,ID];
ID->Association@MapThread[ToString[#1]->#2&,{Columns,#//Rest}]
)&/@Rest[Out];

(*CheckQ[x_]:=Length[x]\[Equal]1;
FFun[x_?CheckQ]:=First[x];

Out=Merge[Out,FFun];
Out=Association@((Normal@Out)/.(HoldPattern[name_\[Rule]FFun[val_]]\[RuleDelayed]Sequence@@MapIndexed[ToString[name]<>"-"<>ToString[First[#2]]\[Rule]#1&,val]));*)
Out=DataTable[Dataset[Association[Out]]];

Out
)
]


(* ::Input::Initialization:: *)
ProcessingQ[f_]:=Or@@(((Head[f]===#||f===#))&/@{Function,Map,MapAt,MapIndexed,KeyValueMap,Apply,ReplacePart,ReplaceAll,Total,Mean,Median,Variance,StandardDeviation});

DTPP[arg_Dataset]:=DataTable[arg];
DTPP[arg_]:=arg;

DTPQ[ds_,ops___]:=Module[
{OPs},
(
OPs=List@ops;
OPs=Replace[OPs,{
s:StringExpression[___]:>KeyFind[s],
Not[s:StringExpression[___]]:>KeyExcept[s]
},{1}];

If[Or@@ProcessingQ/@OPs,
DTPP[DVEvaluate[ds[Sequence@@OPs]]],
DTPP[ds[Sequence@@OPs]]
]
)
];

Dataset[DataTable[ds_Dataset]]^:=ds;

Normal[DataTable[ds_Dataset]]^:=DTPP[Normal[ds]];
DataTable[ds_Dataset][[spec___]]^:=DTPP[ds[[spec]]];

UpValues[DataTable]=Join[UpValues[DataTable],
#/.{
Verbatim[Dataset`Integration`PackagePrivate`ds_Dataset?Dataset`ValidDatasetQ]->DataTable[ds_Dataset],
Dataset`Integration`PackagePrivate`ds->ds,
Dataset`Integration`PackagePrivate`rest->rest,
Dataset`Integration`PackagePrivate`key->key,
Dataset`Integration`PackagePrivate`first->first,
Dataset`Integration`PackagePrivate`arg->arg,
Dataset`EvaluateNormal->DTPP
}&/@Cases[UpValues[Dataset],HoldPattern[((*Except[Verbatim[HoldPattern][((Map|MapAt|MapIndexed|KeyValueMap)[___]|(Map|MapAt|MapIndexed|KeyValueMap)[___][___])]]*)_:>Dataset`EvaluateNormal[_])|(_:>_[Normal[_],___])]]
];

(*Join[DataTable[ds_Dataset],rest__]^:=DTPP[Join[ds,rest]];
Join[left__,DataTable[ds_Dataset],right___]^:=DTPP[Join[left,ds,right]];*)
(*Join[ms:DataTable[ds_Dataset]]^:=ms;*)
JoinAcross[DataTable[ds1_Dataset],DataTable[ds2_Dataset],rest___]^:=DTPP[JoinAcross[ds1,ds2,rest]];

Map[first_,DataTable[ds_Dataset],rest___]^:=DTPQ[Map[first,ds,rest]];
MapAt[first_,DataTable[ds_Dataset],rest___]^:=DTPQ[MapAt[first,ds,rest]];
MapIndexed[first_,DataTable[ds_Dataset],rest___]^:=DTPQ[MapIndexed[first,ds,rest]];
KeyValueMap[first_,DataTable[ds_Dataset],rest___]^:=DTPQ[KeyValueMap[first,ds,rest]];
Apply[first_,DataTable[ds_Dataset],rest___]^:=DTPQ[Apply[first,ds,rest]];
ReplacePart[DataTable[ds_Dataset],rest___]^:=DTPQ[ReplacePart[ds,rest]];
DataTable[ds_Dataset]/.rest___^:=DTPQ[ds/.rest];

Map[arg_][DataTable[ds_Dataset]]^:=DTPQ[Map[arg,ds]];
MapIndexed[arg_][DataTable[ds_Dataset]]^:=DTPQ[MapIndexed[arg,ds]];
KeyValueMap[arg_][DataTable[ds_Dataset]]^:=DTPQ[KeyValueMap[arg,ds]];
Apply[arg_][DataTable[ds_Dataset]]^:=DTPQ[Apply[arg,ds]];
ReplacePart[arg_][DataTable[ds_Dataset]]^:=DTPQ[ReplacePart[arg,ds]];
ReplaceAll[arg_][DataTable[ds_Dataset]]^:=DTPQ[ReplaceAll[arg,ds]];

Total[DataTable[ds_Dataset],rest___]^:=DVEvaluate[Total[Normal[ds],rest]];
Mean[DataTable[ds_Dataset]]^:=DVEvaluate[Mean[Normal[ds]]];
Median[DataTable[ds_Dataset]]^:=DVEvaluate[Median[Normal[ds]]];
Variance[DataTable[ds_Dataset]]^:=DVEvaluate[Variance[Normal[ds]]];
StandardDeviation[DataTable[ds_Dataset]]^:=DVEvaluate[StandardDeviation[Normal[ds]]];

DataTable/:Query[ops___][DataTable[ds_Dataset]]:=DTPQ[ds,ops];
DataTable[ds_Dataset][ops___]:=DTPQ[ds,ops];


(* ::Input::Initialization:: *)
Grid[dt_DataTable,args___]^:=Grid[ToArray@dt,args];
Depth[DataTable[ds_Dataset]]^:=Depth[ds];

Join[dt0_DataTable,Longest[dts__DataTable],n_Integer:1]^:=(
DataTable@Join[Sequence@@(Dataset/@Prepend[List@dts,dt0]),n]
);

JoinT::usage="JoinT[\!\(\*SubscriptBox[\(data\), \(1\)]\),\!\(\*SubscriptBox[\(data\), \(2\)]\),...] returns a DataTable, Dataset, Association or List which results from first transposing the \!\(\*SubscriptBox[\(data\), \(i\)]\), then using Join, and transposing the resulting data object again.";JoinT[dt0_DataTable,Longest[dts__DataTable],n_Integer:1]:=DataTable@Normal@Transpose@Join[Sequence@@(Transpose/@Prepend[List@dts,dt0]),n];
JoinT[dt0_DataTable]:=dt0;

JoinT[dt0_Dataset,Longest[dts__Dataset],n_Integer:1]:=Dataset@Normal@Transpose@Join[Sequence@@(Transpose/@Prepend[List@dts,dt0]),n];

JoinT[dt0_Association,Longest[dts__Association],n_Integer:1]:=Transpose@Join[Sequence@@(Transpose/@Prepend[List@dts,dt0]),n];
JoinT[dt0_List,Longest[dts__List],n_Integer:1]:=Transpose@Join[Sequence@@(Transpose/@Prepend[List@dts,dt0]),n];


(* ::Input::Initialization:: *)
ToArray::usage="ToArray[data] returns a two dimensional List object of the DataTable or Dataset.";
ToArray[ds_Dataset]:=Module[
{Dat,Columns,Lvls,TmpDat},(
Dat=Normal@ds;
Lvls=Depth[DVMeans@Dat];

Which[
Lvls<3,
Prepend[
MapThread[{#1,#2}&,{Keys@Normal@ds,Values@Normal@ds}],
{"ID","Value"}
],

Lvls>=3,(
Columns=DeleteDuplicates@Flatten[Values@Normal@((Keys)/@Dat)];
Prepend[Values@MapIndexed[(Prepend[Values@Normal@#1,Normal[#2][[1,1]]])&,Normal@ds[All,Columns]],Prepend[Columns,"ID"]]
)
]
)];
ToArray[DataTable[ds_Dataset]]:=ToArray[ds];


(* ::Input::Initialization:: *)
DefineItemErrors[data_Association,errors_Association]:=Module[
{OutA},
(
OutA=data;
KeyValueMap[(OutA[#1]=If[NumericQ[OutA[#1]],NDV[OutA[#1],#2],OutA[#1]];)&,errors];
OutA
)
]


(* ::Input::Initialization:: *)
(*If[$VersionNumber\[GreaterEqual]11,
MakeBoxes[tbl_DataTable,StandardForm]^:=(
If[Length[tbl]>19,
(
Column[{
Row[{
Button["Show in Excel",DataTableView@tbl],
Button["Show using \[PlusMinus]",DataTableView[tbl,True]]
},ImageSize\[Rule]{All,Automatic}],
(tbl//First)/.(x_EV\[RuleDelayed]ToString[x])
}]//ToBoxes
),
((tbl//First)/.(x_EV\[RuleDelayed]ToString[x]))//ToBoxes
]
),
MakeBoxes[tbl_DataTable,StandardForm]^:=(
If[Length[tbl]>19,
(
Column[{
Row[{
Button["Show in Excel",DataTableView@tbl],
Button["Show using \[PlusMinus]",DataTableView[tbl,True]]
},ImageSize\[Rule]{All,Automatic}],
(tbl//First)
}]//ToBoxes
),
(tbl//First)//ToBoxes
]
)
]*)

MakeBoxes[tbl_DataTable,StandardForm]^:=(
Column[
{
Row[Join[
If[$UseExcel,
{
Button["Show in Excel",DataTableViewExcel@PrepareListData@tbl],
Button["Show using \[PlusMinus]",DataTableViewExcel@PrepareListData[tbl,True]]
},
{}
],
{
Button["Copy TSV",CopyToClipboard[ExportString[PrepareListData@tbl,"TSV"]]],
Button["Copy using \[PlusMinus]",CopyToClipboard[ExportString[PrepareListData[tbl,True],"TSV"]]]
}
],ImageSize->{All,Automatic}],
(Dataset[tbl])
}
]//ToBoxes

(*If[Length[tbl]>19,
(
Column[{
Row[{
Button["Show in Excel",DataTableView@tbl],
Button["Show using \[PlusMinus]",DataTableView[tbl,True]]
},ImageSize\[Rule]{All,Automatic}],
(tbl//First)
}]//ToBoxes
),
(tbl//First)//ToBoxes
]*)
)


(* ::Input::Initialization:: *)
$ExcelInstance;
$ExcelWorkbook;
$ExcelWorksheet;
$ExcelListObject;
$ExcelCharacters=CharacterRange["A","Z"];
$ExcelColumnRange=Flatten@Prepend[Module[{c1=#},(c1<>#&/@$ExcelCharacters)]&/@$ExcelCharacters,$ExcelCharacters];
$UseExcel=True;

Quiet[
If[!NETObjectQ[$ExcelInstance]||(
$ExcelInstance@Visible;
GetNETException[]!=Null
),
$ExcelInstance=CreateCOMObject["Excel.Application"];
If[!NETObjectQ[$ExcelInstance],
$UseExcel=False;
];
];
];

PrepareListData[tbl_DataTable,separateErrors_:False]:=Module[
{
DataList,Test,HeaderList,i,TempWState
},
(
If[separateErrors,
Module[{ECols},(
DataList=(tbl//ToArray)/.{
(*NDV[v_,e_,_]\[RuleDelayed]{v,e*$ZValue},
DV[v_,e_,__]\[RuleDelayed]{v,e},*)

NDV[v_,se_,n_]:>{v,se*$ZValue},
DV[v_,se_,{rl_?PosNumPatternQ,ru_?PosNumPatternQ},n_]:>{v,{rl,ru}},
DV[v_,se_,i_?PosNumPatternQ,n_]:>{v,i},
DV[v_,se_,n_]:>{v,se*$ZValue}
};
DataList=Map[Switch[#,
_Real,#,
_Integer,#,
_List,#,
_Missing,Null,
_DateObject,DateString[#,{"Date"," ","Time"}],
_,ToString[#]
]&,DataList,{2}];
ECols=Or@@@(Map[ListQ[#]&,DataList,{2}]//Transpose);

DataList[[1]]=Flatten[MapThread[If[#1,{#2,"\[Sigma]("<>#2<>")"},#2]&,{ECols,DataList//First}]];
DataList[[2;;]]=Map[MapThread[If[#1&&!ListQ[#2],{#2,Null},#2]&,{ECols,#}]&,DataList[[2;;]]];
DataList=Flatten/@DataList;
)],
(
DataList=Map[Switch[#,
_Real,#,
_Integer,#,
_Missing,Null,
_DateObject,DateString[#,{"Date"," ","Time"}],
_,ToString[If[AssociationQ[#],EvaluateAssocValues[#/.x_NDV:>ToString[x]],#/.x_NDV:>ToString[x]]]
]&,tbl//ToArray,{2}];
)
];

DataList

(*TempWState=$ExcelWorkbook@Windows@Item[1]@WindowState;
$ExcelWorkbook@Windows@Item[1]@WindowState=-4140;
$ExcelWorkbook@Windows@Item[1]@WindowState=TempWState;*)
)
]


(* ::Input::Initialization:: *)
DataTableViewExcel[data_List]:=Module[
{DataList,HeaderList},
(
DataList=data;
HeaderList=DataList[[1]];
DataList[[1]]=Null&/@DataList[[1]];

Quiet[
If[!NETObjectQ[$ExcelInstance]||(
$ExcelInstance@Visible;
GetNETException[]!=Null
),
$ExcelInstance=CreateCOMObject["Excel.Application"];
If[!NETObjectQ[$ExcelInstance],
Return[$Failed];
];
];
$ExcelInstance@Visible=True;

If[!NETObjectQ[$ExcelWorkbook]||$ExcelWorkbook@Worksheets==$Failed||(
$ExcelWorkbook@Worksheets;
GetNETException[]!=Null
),
$ExcelWorkbook = $ExcelInstance@Workbooks@Add[];
];

If[NETObjectQ[$ExcelWorksheet]&&($ExcelWorksheet@Visible==$Failed||(
$ExcelWorksheet@Visible;
GetNETException[]!=Null
)),
ReleaseNETObject[$ExcelWorksheet];
$ExcelWorksheet = $ExcelWorkbook@Worksheets@Add[];
];

If[!NETObjectQ[$ExcelWorksheet],
ReleaseNETObject[$ExcelWorksheet];
$ExcelWorksheet = $ExcelWorkbook@Worksheets@Item[1];
$ExcelInstance@Visible = True;

$ExcelWorksheet@Range["B2"]@Select[];
$ExcelWorkbook@Windows@Item[1]@FreezePanes=True;
];

If[NETObjectQ[$ExcelListObject],$ExcelListObject@Delete[]];
ReleaseNETObject[$ExcelListObject];

];

$ExcelWorksheet@Range["B2","ZZ2"]@EntireColumn@HorizontalAlignment=4;
$ExcelWorksheet@Range["A1","A2"]@EntireColumn@HorizontalAlignment=1;
$ExcelWorksheet@Range["A1","ZZ1"]@EntireRow@HorizontalAlignment=3;

$ExcelWorksheet@Range["A1",$ExcelColumnRange[[DataList//First//Length]]<>ToString[DataList//Length]]@Value=DataList;

$ExcelWorksheet@Range["A1",$ExcelColumnRange[[DataList//First//Length]]<>"1"]@Value=HeaderList;

$ExcelWorksheet@Range["A1",$ExcelColumnRange[[DataList//First//Length]]<>ToString[DataList//Length]]@Columns@AutoFit[];

$ExcelListObject=$ExcelWorksheet@ListObjects@Add[1,$ExcelWorksheet@Range["A1",$ExcelColumnRange[[DataList//First//Length]]<>ToString[DataList//Length]],False,1];
$ExcelListObject@TableStyle=$ExcelWorkbook@TableStyles@Item[23];
$ExcelListObject@ShowAutoFilter=False;
)
]


(* ::Input::Initialization:: *)
KeyFind::usage="KeyFind[str] represents an operator which can be used to select elements whose key is the provided string.
KeyFind[patt] represents an operator which can be used to select elements whose key is the provided string pattern.";
KeyFind[str_]:=KeySelect[StringMatchQ[ToString[#],str]&];
KeyFind[str_String]:=Key[str];
KeyFind[assoc_Association,str_]:=KeySelect[assoc,StringMatchQ[ToString[#],str]&];

KeyExcept::usage="KeyExcept[str] represents an operator which can be used to select elements whose key is not the provided string.
KeyExcept[patt] represents an operator which can be used to select elements whose key is not matched by the provided string pattern.";
KeyExcept[str_]:=KeySelect[StringFreeQ[ToString[#],str]&];
KeyExcept[str_String]:=KeyDrop[str];
KeyExcept[assoc_Association,str_]:=KeySelect[assoc,StringFreeQ[ToString[#],str]&];

Unprotect[StringExpression];
DataTable[ds_][op1_StringExpression,op2_StringExpression,opts:OptionsPattern[]]^:=DataTable[ds][KeyFind[op1],KeyFind[op2],opts];
DataTable[ds_][op1_StringExpression,op2_,opts:OptionsPattern[]]^:=DataTable[ds][KeyFind[op1],op2,opts];
DataTable[ds_][op1_,op2_StringExpression,opts:OptionsPattern[]]^:=DataTable[ds][op1,KeyFind[op2],opts];
DataTable[ds_][op1_StringExpression,opts:OptionsPattern[]]^:=DataTable[ds][KeyFind[op1],opts];
Protect[StringExpression];


(* ::Input::Initialization:: *)
Unprotect[Not,StringExpression];
DataTable[ds_][op1_Not,op2_StringExpression,opts:OptionsPattern[]]^:=DataTable[ds][KeyExcept[op1//First],KeyFind[op2],opts];
DataTable[ds_][op1_Not,op2_,opts:OptionsPattern[]]^:=DataTable[ds][KeyExcept[op1//First],op2,opts];
DataTable[ds_][op1_Not,opts:OptionsPattern[]]^:=DataTable[ds][KeyExcept[op1//First],opts];
Protect[Not,StringExpression];


(* ::Input::Initialization:: *)
NV::usage="NV[datatable] returns the data (not the keys) in datatable, converted to a normal expression.";
NV[ds_Dataset]:=Normal@Values@Map[If[AssociationQ[#],Normal@Values[#],#]&,ds,Infinity];
NV[DataTable[ds_Dataset]]:=NV[ds];

FV::usage="FV[datatable] returns the first value in datatable.";
FV[ds_Dataset]:=First@Map[If[AssociationQ[#]||ListQ[#],First[#],#]&,ds,Infinity];
FV[DataTable[ds_Dataset]]:=FV[ds];

ColNames::usage="ColNames[datatable] returns the list of the columns in the data table.";
ColNames[ds_Dataset]:=If[Depth[DVMeans@Normal@ds]<3,{},DeleteDuplicates@Flatten[Values@Normal@(Keys/@Normal@ds)]];
ColNames[DataTable[ds_Dataset]]:=ColNames[ds];

NK::usage="NK[datatable] returns the list of the row keys (IDs) in the data table."
NK[ds_Dataset]:=Normal@Keys@ds;
NK[DataTable[ds_Dataset]]:=NK[ds];


(* ::Input::Initialization:: *)
MeanAndStandardError::usage="MeanAndStandardError represents an an operator which can be applied to DataTable objects to calculate the means and standard errors, represented by an NDV object, e.g. of all rows.";
MeanAndStandardError=Merge[Values@(#),NDV[Mean[#],StandardDeviation[#]/Sqrt[Length[#]]//N]&]&;

MeanAndStandardDeviation::usage="MeanAndStandardDeviation represents an an operator which can be applied to DataTable objects to calculate the means and standard deviations, represented by an NDV object, e.g. of all rows.";
MeanAndStandardDeviation=Merge[Values@(#),NDV[Mean[#],StandardDeviation[#]//N]&]&;


(* ::Input::Initialization:: *)
FixedEffectMetaAnalysis::usage="FixedEffectMetaAnalysis represents an an operator which can be applied to DataTable objects to carry out a fixed effect meta-analysis.";
FixedEffectMetaAnalysis=Merge[Values@#,Module[{WDat},
(
WDat=WeightedData[DVMeans@#,1/#^2&/@DVSEs[#]];
NDV[Mean[WDat],Sqrt[1/(Total[1/#^2&/@DVSEs[#]])]//N]
)
]&]&;


(* ::Input::Initialization:: *)
RandomEffectsMetaAnalysis::usage="RandomEffectsMetaAnalysis represents an an operator which can be applied to DataTable objects to carry out a random effect meta-analysis.";
RandomEffectsMetaAnalysis=Merge[Values@#,Module[{WT,WTS,ESWT,ESSWT,Cv,Q,TS,Wtot,ESWtot},
(
WT=(1/(DVSEs[#]^2)^2)&/@#;
WTS=(1/(DVSEs[#]^2)^2)^2&/@#;
ESWT=(DVMeans[#]/(DVSEs[#]^2)^2)&/@#;
ESSWT=(DVMeans[#]^2/(DVSEs[#]^2)^2)&/@#;
Q=Total[ESSWT]-Total[ESWT]^2/Total[WT];
Cv=Total[WT]-Total[WTS]/Total[WT];

TS=If[Q>Length[#]-1,(Q-(Length[#]-1))/Cv,0];

Wtot=1/((DVSEs[#]^2)+TS)&/@#;
ESWtot=DVMeans[#]/((DVSEs[#]^2)+TS)&/@#;

NDV[Total[ESWtot]/Total[Wtot],Sqrt[1/Total[Wtot]]//N]
)
]&]&;


(* ::Input::Initialization:: *)
GroupByKey::usage="GroupByKey[datatable,f] returns a DataTable in which the rows in in the input data table have been grouped into sublists with distinct f[row key (ID)].";
GroupByKey[data_DataTable,f_]:=Module[
{Dat},
(
Dat=data[MapIndexed[Prepend[#1,"GroupingKey"->f[#2[[1,1]]]]&]];
Dat[GroupBy[#,"GroupingKey"/.#&]&][All,All,KeyDrop["GroupingKey"]]
)
];
GroupByKey[data_DataTable,f_,2]:=Module[
{Dat},
(
Dat=data[All,MapIndexed[Prepend[{"Value"->#1},"GroupingKey"->f[#2[[1,1]]]]&,#,1]&];
Dat[All,GroupBy[#,"GroupingKey"/.#&]&][All,All,All,"Value"/.#&]
)
]


(* ::Input::Initialization:: *)
PartitionAssociation::usage="PartitionAssociation[assoc,n] partitions the association assoc into nonoverlapping subassociations of length n and returns a new association.
PartitionAssociation[assoc,n,d] generates subassociations with offet d.
PartitionAssociation[assoc,n,d,i] generates subassociations with offet d and the keys of the resulting association are generated from item number i of the partition group.
";
PartitionAssociation[assoc_,n_,d_,i_:1]:=Association[(#[[i,1]]->Association@#)&/@Partition[Normal@assoc,n,d]]
PartitionAssociation[assoc_,n_]:=Association[(#[[1,1]]->Association@#)&/@Partition[Normal@assoc,n]]


(* ::Input::Initialization:: *)
SetAttributes[QueryTo,HoldFirst];
Options[QueryTo]=Options[Query];
QueryTo::usage="QueryTo[obj,\!\(\*SubscriptBox[\(operator\), \(1\)]\),\!\(\*SubscriptBox[\(operator\), \(2\)]\),...] applies the query with the given operators to the Dataset or Datatable obj.";
QueryTo[dt_,op___,opts:OptionsPattern[]]:=dt=Query[op,opts][dt]


(* ::Input::Initialization:: *)
ReplaceDatum::usage="ReplaceDatum[obj,row,column,value] replaces the value of the datum in the row with a certain key (ID) and the column with a certain key with a new value.";
ReplaceDatum[dt_,row_,column_,val_]:=(
dt[{row->(<|#,column->val|>&)}]
)


(* ::Input::Initialization:: *)
FoldValues::usage="FoldValues[f,x,assoc] applies FoldList[f,x,list] to the values of items in an association and returns a new association using the keys from the second item in each operation. The arguments x is optional.
FoldValues[f] represents the operator form of FoldValues.";
FoldValues[f_,x_,assoc_Association]:=<|FoldList[#2[[1]]->f[#1[[2]],#2[[2]]]&,If[Head[x[First@Values@assoc]]===Rule,x[First@Values@assoc],"Key0"->x[First@Values@assoc]],Normal@assoc]|>;
FoldValues[f_,assoc_Association]:=<|FoldList[#2[[1]]->f[#1[[2]],#2[[2]]]&,Normal@assoc]|>;
FoldValues[f_]:=<|FoldList[#2[[1]]->f[#1[[2]],#2[[2]]]&,Normal@#]|>&;


(* ::Input::Initialization:: *)
$DefaultPlotColors=First/@("DefaultPlotStyle"/.(Method/.Charting`ResolvePlotTheme[Automatic,ListPlot]));

$DefaultPlotMarkerG={
{EdgeForm[Black],Disk[]},
{EdgeForm[Black],Rectangle[]},
{EdgeForm[Black],Polygon[{{1,0},{0,Sqrt[3]},{-1,0}}]},
{EdgeForm[Black],Polygon[{{0.,1.},{-1.,0.},{0.,-1.},{1.,0.}}]},
{EdgeForm[Black],Polygon[{{1,0},{0,-Sqrt[3]},{-1,0}}]},
{EdgeForm[Black],Polygon[{{0,0},{-1,1},{0,0},{-1,-1},{0,0},{1,-1},{0,0},{1,1},{0,0}}]},
{EdgeForm[Black],Polygon[{{1,0},{0,Sqrt[3]},{-1,0}}]},
{EdgeForm[Black],Polygon[{{0,0},{Sqrt[3],-1},{Sqrt[3],1}}]},
{EdgeForm[Black],Polygon[{{Sqrt[3],0},{0,-1},{0,1}}]}
};

$DefaultPlotMarkers=(Graphics[#,ImageSize->10]&/@MapIndexed[
Prepend[$DefaultPlotMarkerG[[Mod[#2[[1]],Length[$DefaultPlotMarkerG],1]]],#1]&,$DefaultPlotColors]);

Options[DTPlot]=Options[ErrorListPlot];
SetOptions[DTPlot,PlotMarkers->Automatic];

DTPlot::usage="DTPlot[datatable] returns a plots of the DataTable object.
DTPlot[{\!\(\*SubscriptBox[\(datatable\), \(1\)]\),\!\(\*SubscriptBox[\(datatable\), \(2\)]\),...}] returns a plots of the DataTable objects.";
DTPlot[data_,opts:OptionsPattern[]]:=Module[
{DTs},
(
DTs=If[Head[data]===DataTable,{data},data];
DTs=If[Head[#]===DataTable,DeleteMissing[#,1,1],#]&/@DTs;

ErrorListPlot[
If[Head[#]===DataTable,PrepareErrorList[#,!(OptionValue[PlotMarkers]===None)]&@#,#]&/@DTs,
Join[{opts},{
Frame->True,
FrameStyle->Black,
FrameLabel->(ColNames@First@DTs)[[;;2]],
PlotLegends->Automatic,
PlotRangePadding->Scaled[0.1],
PlotRange->Full,
ImageSize->400,
(*PlotMarkers\[Rule]$DefaultPlotMarkers,
Method\[Rule]{"OptimizePlotMarkers"\[Rule]False},*)
PlotStyle->$DefaultPlotColors,
FrameTicksStyle->12,
LabelStyle->12,
Axes->False,
ErrorBarFunction->(ErrorBarPlots`Private`ebarfun[##]/.l_Line:>{Black,Thin,Dashing[{}],l}&)
}]]/.GraphicsComplex[x_,y_,z___]:>(GraphicsComplex[x,Replace[y,k_List:>If[(k//First)===Opacity[0.`],k,Reverse@k],{3}],z])
)
]


(* ::Input::Initialization:: *)
PrepareErrorList[dt_DataTable,showErrors_]:=Module[
{},
(
Prepend[MapThread[Tooltip[{
DVMeans@{#1[[1]],#1[[2]]},
ErrorBar[
If[showErrors,Switch[Head[#1[[1]]],
NDV,IV[#1[[1]]][[2]],
IV,#1[[1,2]],
DV,If[#1[[1,3]]===Null,0,#1[[1,3]]],
_,0
],0],
If[showErrors,If[showErrors,Switch[Head[#1[[2]]],
NDV,IV[#1[[2]]][[2]],
IV,#1[[2,2]],
DV,If[#1[[2,3]]===Null,0,#1[[2,3]]],
_,0
],0],0]
]
},ToString[#2]<>" "<>ToString[ToString/@#1]]&,
{NV@dt,NK@dt}
],{Null,Null}]
)
]


(* ::Input::Initialization:: *)
DefineStyle["DataAnalysis",
{
FigObject->{FontFamily->"Helvetica",FontSize->12},
Figure->{CanvasMargin->{{0.8,0.3},{0.5,0.5}}},
FigurePanel->{TickFontSize->11,FrameFontSize->12,PanelLetterFontSize->12,ExtendRange->0.07},
Multipanel->{TickFontSize->11,FrameFontSize->12,PanelLetterFontSize->12,ExtendRange->0.07},
DataLegend->{FontSize->12},
FigLabel->{FontSize->12}
}
];

$SciDrawPlotMarkers={"Circle","Square","Diamond","UpTriangle","DownTriangle","LeftTriangle","RightTriangle"};
MapIndexed[(
DefineStyle[{"DataAnalysisCurves",#2[[1]]},{
DataLine->{
LineColor->#1,
LineThickness->Medium,
LineDashing->None
},
FigLine->{
LineColor->#1,
LineThickness->Medium,
LineDashing->None
},
DataSymbol->{
SymbolShape->$SciDrawPlotMarkers[[Mod[#2[[1]],Length[$SciDrawPlotMarkers]]/.{0->Length[$SciDrawPlotMarkers]}]],FillColor->#1,SymbolSize->8
}
}];
)&,$DefaultPlotColors];


(* ::Input::Initialization:: *)
Options[SciDrawPlot]={XErrorBars->True,YErrorBars->True(*,XLogScale\[Rule]False,YLogScale\[Rule]False*),AnchorName->None}~Join~Options[DataPlot];
XErrorBars::usage="XErrorBars is an option that specifies if error bars should be drawn for the x-axis.";
YErrorBars::usage="YErrorBars is an option that specifies if error bars should be drawn for the y-axis.";

SciDrawPlot::usage="SciDrawPlot[{{\!\(\*SubscriptBox[\(x\), \(1\)]\),\!\(\*SubscriptBox[\(y\), \(1\)]\)},{\!\(\*SubscriptBox[\(x\), \(2\)]\),\!\(\*SubscriptBox[\(y\), \(2\)]\)},...}] returns a DataPlot object which can be used to plot the data inside a SciDraw Figure. Error bars can be plotted if \!\(\*SubscriptBox[\(x\), \(i\)]\) and \!\(\*SubscriptBox[\(y\), \(i\)]\) are NDV-, DV-, or IV-type objects.
SciDrawPlot[datatable] returns a DataPlot object with data in a DataTable, which can be used inside a SciDraw Figure.";
SciDrawPlot[dt_List,opts:OptionsPattern[]]:=Module[
{Data,XAsymErrors=False,YAsymErrors=False},
(
Data=dt;

If[OptionValue[XErrorBars],
(
XAsymErrors=TrueQ[Or@@(MatchQ[#,DV[_,_,{_?PosNumPatternQ,_?PosNumPatternQ},_]|IV[_,{_?PosNumPatternQ,_?PosNumPatternQ}]]&/@Data[[All,1]])];
Data=If[XAsymErrors,
Switch[Head[#[[1]]],
NDV,{#[[1,1]],#[[2]],IV[#[[1]]][[2]],IV[#[[1]]][[2]]},
DV,{#[[1,1]],#[[2]],If[#[[1,3]]===Null,{0,0},If[ListQ[#[[1,3]]],#[[1,3]],{#[[1,3]],#[[1,3]]}]]},
IV,{#[[1,1]],#[[2]],If[ListQ[#[[1,2]]],#[[1,2]],{#[[1,2]],#[[1,2]]}]},
_,{#[[1]],#[[2]],0,0}
]&/@Data,
Switch[Head[#[[1]]],
NDV,{#[[1,1]],#[[2]],IV[#[[1]]][[2]]},
DV,{#[[1,1]],#[[2]],If[#[[1,3]]===Null,0,#[[1,3]]]},
IV,{#[[1,1]],#[[2]],#[[1,2]]},
_,{#[[1]],#[[2]],0}
]&/@Data
];
Data=Flatten/@Data;
),
Data[[All,1]]=Data[[All,1]]//DVMeans;
];

If[OptionValue[YErrorBars],
(
YAsymErrors=TrueQ[Or@@(MatchQ[#,DV[_,_,{_?PosNumPatternQ,_?PosNumPatternQ},_]|IV[_,{_?PosNumPatternQ,_?PosNumPatternQ}]]&/@Data[[All,2]])];

Data=If[YAsymErrors,
Switch[Head[#[[2]]],
NDV,If[Length[#]>2,
{#[[1]],#[[2,1]],#[[3;;]],IV[#[[2]]][[2]],IV[#[[2]]][[2]]},
{#[[1]],#[[2,1]],IV[#[[2]]][[2]],IV[#[[2]]][[2]]}],
DV,If[Length[#]>2,
{#[[1]],#[[2,1]],#[[3;;]],If[#[[2,3]]===Null,{0,0},If[ListQ[#[[2,3]]],#[[2,3]],{#[[2,3]],#[[2,3]]}]]},
{#[[1]],#[[2,1]],If[#[[2,3]]===Null,{0,0},If[ListQ[#[[2,3]]],#[[2,3]],{#[[2,3]],#[[2,3]]}]]}
],
IV,If[Length[#]>2,
{#[[1]],#[[2,1]],#[[3;;]],If[ListQ[#[[2,2]]],#[[2,2]],{#[[2,2]],#[[2,2]]}]},
{#[[1]],#[[2,1]],If[ListQ[#[[2,2]]],#[[2,2]],{#[[2,2]],#[[2,2]]}]}
],
_,If[Length[#]>2,
{#[[1]],#[[2]],#[[3;;]],0,0},
{#[[1]],#[[2]],0,0}
]
]&/@Data,
Switch[Head[#[[2]]],
NDV,If[Length[#]>2,
{#[[1]],#[[2,1]],#[[3;;]],IV[#[[2]]][[2]]},
{#[[1]],#[[2,1]],IV[#[[2]]][[2]]}],
DV,If[Length[#]>2,
{#[[1]],#[[2,1]],#[[3;;]],If[#[[2,3]]===Null,0,#[[2,3]]]},
{#[[1]],#[[2,1]],If[#[[2,3]]===Null,0,#[[2,3]]]}
],
IV,If[Length[#]>2,
{#[[1]],#[[2,1]],#[[3;;]],#[[2,2]]},
{#[[1]],#[[2,1]],#[[2,2]]}
],
_,If[Length[#]>2,
{#[[1]],#[[2]],#[[3;;]],0},
{#[[1]],#[[2]],0}
]
]&/@Data
];
Data=Flatten/@Data;
),
Data[[All,2]]=Data[[All,2]]//DVMeans;
];

(*Data=MapThread[Tooltip[#1,#2<>" "<>ToString[{
ToString@If[TrueQ[OptionValue[XLogScale]],If[#1[[3]]>0,EVE@(10^EV[#1[[1]],#1[[3]]]),10^#1[[1]]],#1[[1]]],
ToString@If[TrueQ[OptionValue[YLogScale]],If[TrueQ[XErrorBars],
If[#1[[3]]>0,EVE@(10^EV[#1[[1]],#1[[3]]]),10^#1[[1]]],
If[#1[[4]]>0,EVE@(10^EV[#1[[1]],#1[[4]]]),10^#1[[1]]]
],#1[[2]]]
}]]&,{Data,
labels}];*)

If[TrueQ[OptionValue[AnchorName]!=None],
DataPlot[[OptionValue[AnchorName]]][Data,Sequence@@DeleteCases[Join[{
If[OptionValue[XErrorBars],XErrorColumn->If[XAsymErrors,{3,4},3],Null],
If[OptionValue[YErrorBars],YErrorColumn->If[OptionValue[XErrorBars],If[XAsymErrors,If[YAsymErrors,{5,6},5],If[YAsymErrors,{4,5},4]],If[YAsymErrors,{3,4},3]],Null]
},FilterRules[{opts}, Options[DataPlot]]],Null]],
DataPlot[Data,Sequence@@DeleteCases[Join[{
If[OptionValue[XErrorBars],XErrorColumn->If[XAsymErrors,{3,4},3],Null],
If[OptionValue[YErrorBars],YErrorColumn->If[OptionValue[XErrorBars],If[XAsymErrors,If[YAsymErrors,{5,6},5],If[YAsymErrors,{4,5},4]],If[YAsymErrors,{3,4},3]],Null]
},FilterRules[{opts}, Options[DataPlot]]],Null]]
]
)
];

SciDrawPlot[dt_DataTable,opts:OptionsPattern[]]:=Module[
{TrimmedData},
(
TrimmedData=DeleteMissing[dt,1,1];
SciDrawPlot[TrimmedData//NV,(*TrimmedData//NK,*)Join[{AnchorName->ColNames[dt][[2]]},{opts}]]
)];


(* ::Input::Initialization:: *)
SetAttributes[Fig,HoldFirst];
Options[Fig]=Join[{XLogScale->False,YLogScale->False,Joined->False,PlotMarkers->True,ShowLegend->Automatic,PlotLegends->Automatic},Options[FigurePanel],Options[Figure],Options[SciDrawPlot]];
XLogScale::usage="XLogScale is an option that specifies if Log10[x] should be applied to the x coordinates in the plot.";
YLogScale::usage="YLogScale is an option that specifies if Log10[y] should be applied to the y coordinates in the plot.";
ShowLegend::usage="ShowLegend is an option that specifies if a plot legend should be shown.";

Fig::usage="Fig[{\!\(\*SubscriptBox[\(object\), \(1\)]\),\!\(\*SubscriptBox[\(object\), \(2\)]\),...}] returns a figure with plots of the objects using the SciDraw package. The objects can be of the type DataTable or Graphic.";
Fig::nodatatable="At least one DataTable must be displayed if no PlotRange option is not specified.";
Fig[plots_,opts:OptionsPattern[]]:=Module[
{TmpStyles,DTs,PR,ColLabels},
(
DTs=If[Head[plots]===DataTable,{plots},Cases[plots,_DataTable]];
DTs=DeleteMissing[#,1,1]&/@DTs;
PR=If[ListQ[OptionValue[PlotRange]],OptionValue[PlotRange],CalcPlotRange[DTs,FilterRules[{opts}, Options[CalcPlotRange]]]];
ColLabels=ColNames@First@DTs;

If[TrueQ[OptionValue[XLogScale]]||TrueQ[OptionValue[YLogScale]],
DTs=#[All,<|
"x"->If[TrueQ[OptionValue[XLogScale]],Log10@#[[1]],#[[1]]],
"y"->If[TrueQ[OptionValue[YLogScale]],Log10@#[[2]],#[[2]]]
|>&]&/@DTs;
];

Figure[
FigurePanel[
{
SetOptions[DataPlot,DataSymbol->{Show->OptionValue[PlotMarkers]}],
SetOptions[DataPlot,DataLine->{Show->OptionValue[Joined]}],

(FigGraphics/@Cases[plots,_Graphics]),
(ReleaseHold/@Cases[plots,Except[_DataTable|_Graphics]]),

MapIndexed[SciDrawPlot[#1,Append[FilterRules[{opts}, Options[SciDrawPlot]],
Style->{"DataAnalysisCurves",Mod[#2[[1]],Length[$DefaultPlotColors],1]}]]&,DTs](*,*)

If[(OptionValue[ShowLegend]===Automatic&&Length[DTs]>1)||ListQ[OptionValue[ShowLegend]],
DataLegend[If[OptionValue[ShowLegend]===Automatic,Scaled[{0.8,0.2}],Scaled[OptionValue[ShowLegend]]],
If[OptionValue[PlotLegends]===Automatic,
Table[{{"DataAnalysisCurves",Mod[i,Length[$DefaultPlotColors],1]},ToString[i]},{i,1,Length[DTs]}],
Table[{{"DataAnalysisCurves",Mod[i,Length[$DefaultPlotColors],1]},OptionValue[PlotLegends][[i]]},{i,1,Length[DTs]}]
],
FontFamily->"sans-serif"
],
Null
]
},
FilterRules[{opts}, Options[FigurePanel]],
PlotRange->PR,
Frame->True,
XFrameLabel->If[OptionValue[XLogScale],Log[ColLabels[[1]]],ColLabels[[1]]],
YFrameLabel->If[OptionValue[YLogScale],Log[ColLabels[[2]]],ColLabels[[2]]]
],
FilterRules[{opts}, Options[Figure]],
CanvasSize->1*{6,4},
Style->"DataAnalysis"
]
)
]

Options[CalcPlotRange]={XLogScale->False,YLogScale->False(*,PlotRangePadding\[Rule]0.05*)};
CalcPlotRange[dt_DataTable,opts:OptionsPattern[]]:=Module[
{XMin,XMax,YMin,YMax,data},
(
data=DeleteMissing[dt,1,1];

XMin=Min@data[All,1,DVMeans];
If[TrueQ[OptionValue[XLogScale]],XMin=Log10[XMin]];
XMax=Max@data[All,1,DVMeans];
If[TrueQ[OptionValue[XLogScale]],XMax=Log10[XMax]];

(*XMin=XMin-Abs[XMax-XMin]*If[ListQ[OptionValue[PlotRangePadding]],OptionValue[PlotRangePadding][[1]],OptionValue[PlotRangePadding]];
If[XMin\[NotEqual]0,XMin=Floor[XMin,10^((Log10[Abs[XMin]]//Floor)-1)]];
XMax=XMax+Abs[XMax-XMin]*If[ListQ[OptionValue[PlotRangePadding]],OptionValue[PlotRangePadding][[1]],OptionValue[PlotRangePadding]];
If[XMax\[NotEqual]0,XMax=Ceiling[XMax,10^((Log10[Abs[XMax]]//Floor)-1)]];*)

YMin=Min@data[All,2,DVMeans];
If[TrueQ[OptionValue[YLogScale]],YMin=Log10[YMin]];
YMax=Max@data[All,2,DVMeans];
If[TrueQ[OptionValue[YLogScale]],YMax=Log10[YMax]];

(*YMin=YMin-Abs[YMax-YMin]*If[ListQ[OptionValue[PlotRangePadding]],OptionValue[PlotRangePadding][[2]],OptionValue[PlotRangePadding]];
If[YMin\[NotEqual]0,YMin=Floor[YMin,10^((Log10[Abs[YMin]]//Floor)-1)]];
YMax=YMax+Abs[YMax-YMin]*If[ListQ[OptionValue[PlotRangePadding]],OptionValue[PlotRangePadding][[2]],OptionValue[PlotRangePadding]];
If[YMax\[NotEqual]0,YMax=Ceiling[YMax,10^((Log10[Abs[YMax]]//Floor)-1)]];*)

(*If[TrueQ[OptionValue[XLogScale]],XMin=Floor[#,0.1]&@Log10[XMin];XMax=Ceiling[#,0.1]&@Log10[XMax];];
If[TrueQ[OptionValue[YLogScale]],YMin=Floor@Log10[YMin];YMax=Ceiling@Log10[YMax];];*)

{{XMin,XMax},{YMin,YMax}}
)
];

CalcPlotRange[data_List,opts:OptionsPattern[]]:=Module[{Ranges},
(
Ranges=(CalcPlotRange[#,opts]&/@data);
If[Length[Ranges]==0,ErrorMessage[Fig::nodatatable]];
{
{Min@@Ranges[[All,1,1]],Max@@Ranges[[All,1,2]]},
{Min@@Ranges[[All,2,1]],Max@@Ranges[[All,2,2]]}
}
)
];


(* ::Input::Initialization:: *)
FigColumn::usage="FigColumn[datatable,x] returns a figure with vertically arranged panels plotting variable x vs. all other columns in the DataTable using the SciDraw package.
FigColumn[datatable,x,{\!\(\*SubscriptBox[\(ycol\), \(1\)]\),\!\(\*SubscriptBox[\(ycol\), \(2\)]\),..}] returns a figure with vertically arranged panels plotting variable x vs. columns \!\(\*SubscriptBox[\(ycol\), \(i\)]\) in the DataTable using the SciDraw package.
FigColumn[{\!\(\*SubscriptBox[\(dt\), \(1\)]\),\!\(\*SubscriptBox[\(dt\), \(2\)]\),..},x] plots the variables for DataTables \!\(\*SubscriptBox[\(dt\), \(i\)]\).";

Options[FigColumn]=Join[{XLogScale->False,YLogScale->False,Joined->False,PlotMarkers->True,ShowLegend->Automatic,PlotLegends->Automatic},Options[Multipanel],Options[Figure],Options[SciDrawPlot]];

FigColumn[data_,x_,opts:OptionsPattern[]]:=FigColumn[data,x,All,opts];
FigColumn[data_,x_,y_:All,opts:OptionsPattern[]]:=Module[
{YList},
(
YList=If[y===All,If[ListQ[data],Rest@ColNames@data[[1]],Rest@ColNames@data],y];

Figure[
Multipanel[
MapIndexed[(FigurePanel[
{
SetOptions[DataPlot,DataSymbol->{Show->OptionValue[PlotMarkers]}],
SetOptions[DataPlot,DataLine->{Show->OptionValue[Joined]}],

If[ListQ[data],
Module[{YVar=#1},MapIndexed[SciDrawPlot[#1[All,{x,YVar}],Style->{"DataAnalysisCurves",Mod[#2[[1]],Length[$DefaultPlotColors],1]},FilterRules[{opts}, Options[SciDrawPlot]]]&,data]],
{SciDrawPlot[data[All,{x,#1}],Style->{"DataAnalysisCurves",1},FilterRules[{opts}, Options[SciDrawPlot]]]}
],
If[#2[[1]]==1,
If[(OptionValue[ShowLegend]===Automatic&&ListQ[data]&&Length[data]>1)||ListQ[OptionValue[ShowLegend]],
DataLegend[If[OptionValue[ShowLegend]===Automatic,Scaled[{0.8,0.2}],Scaled[OptionValue[ShowLegend]]],
If[OptionValue[PlotLegends]===Automatic,
Table[{{"DataAnalysisCurves",Mod[i,Length[$DefaultPlotColors],1]},ToString[i]},{i,1,If[ListQ[data],Length[data],1]}],
Table[{{"DataAnalysisCurves",Mod[i,Length[$DefaultPlotColors],1]},OptionValue[PlotLegends][[i]]},{i,1,If[ListQ[data],Length[data],1]}]
],
FontFamily->"sans-serif"
],
Null
]
]
},
{#2[[1]],1}
])&,YList],
Dimensions->{Length[YList],1},

FilterRules[{opts}, Options[Multipanel]],
XPlotRange->CalcPlotRange[If[ListQ[data],Flatten[Module[{YVar=#},#[All,{x,YVar}]&/@data]&/@YList],data[All,{x,#}]&/@YList]][[1]],
YPlotRange->If[ListQ[data],CalcPlotRange[Module[{YVar=#},#[All,{x,YVar}]&/@data]][[2]]&/@YList,CalcPlotRange[data[All,{x,#}]][[2]]&/@YList],
Frame->True,
XFrameLabel->x,
YFrameLabel->YList,

PanelLetter->None,
YPanelGaps->0.1
],

FilterRules[{opts}, Options[Figure]],
CanvasSize->1*{6,2*Length[YList]},
Style->"DataAnalysis"
]
)
]


(* ::Input::Initialization:: *)
Options[FigGrid]=Join[{XLogScale->False,YLogScale->False,Joined->False,PlotMarkers->True,ShowLegend->Automatic,PlotLegends->Automatic,FigPanelExtras-><||>},Options[Multipanel],Options[Figure],Options[SciDrawPlot]];

FigPanelExtras::usage="FigPanelExtras is an option for FigGrid which specifies an association with \"extra\" graphic elements to plot in individual grid cells.";

FigGrid::usage="FigGrid[\!\(\*
StyleBox[\"data\",\nFontSlant->\"Italic\"]\),\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)] returns a figure grid with plots of data in DataTable objects, where column \!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\) is plotted on the x-axis, using the SciDraw package.
FigGrid[\!\(\*
StyleBox[\"data\",\nFontSlant->\"Italic\"]\),\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"{\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[SubscriptBox[\"ycol\", \"1\"],\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[SubscriptBox[\"ycol\", \"2\"],\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"..\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"}\",\nFontSlant->\"Italic\"]\)] plots only columns specified by {\!\(\*SubscriptBox[\(ycol\), \(1\)]\),\!\(\*SubscriptBox[\(ycol\), \(2\)]\),..} in the figure grid.";

FigGrid[data_,x_,opts:OptionsPattern[]]:=FigGrid[data,x,All,opts];
FigGrid[data_,x_,y_:All,opts:OptionsPattern[]]:=Module[
{YList},
(
YList=If[y===All,If[ListQ[data[[1]]],DeleteCases[ColNames@data[[1,1]],x],DeleteCases[ColNames@data[[1]],x]],y];

Figure[
Multipanel[
MapIndexed[
Module[{ColData=#1,ColI=#2[[1]]},
(
MapIndexed[(FigurePanel[
{
If[KeyExistsQ[OptionValue[FigPanelExtras],{#1,ColI}],OptionValue[FigPanelExtras][{#1,ColI}]//ReleaseHold],

SetOptions[DataPlot,DataSymbol->{Show->OptionValue[PlotMarkers]}],
SetOptions[DataPlot,DataLine->{Show->OptionValue[Joined]}],

If[ListQ[ColData],
Module[{YVar=#1},MapIndexed[SciDrawPlot[#1[All,{x,YVar}],Style->{"DataAnalysisCurves",Mod[#2[[1]],Length[$DefaultPlotColors],1]},FilterRules[{opts}, Options[SciDrawPlot]]]&,ColData]],
{SciDrawPlot[ColData[All,{x,#1}],Style->{"DataAnalysisCurves",1},FilterRules[{opts}, Options[SciDrawPlot]]]}
],

If[#2[[1]]==1&&ColI==1,
If[(OptionValue[ShowLegend]===Automatic&&ListQ[ColData]&&Length[ColData]>1)||ListQ[OptionValue[ShowLegend]],
DataLegend[If[OptionValue[ShowLegend]===Automatic,Scaled[{0.8,0.2}],Scaled[OptionValue[ShowLegend]]],
If[OptionValue[PlotLegends]===Automatic,
Table[{{"DataAnalysisCurves",Mod[i,Length[$DefaultPlotColors],1]},ToString[i]},{i,1,If[ListQ[ColData],Length[ColData],1]}],
Table[{{"DataAnalysisCurves",Mod[i,Length[$DefaultPlotColors],1]},OptionValue[PlotLegends][[i]]},{i,1,If[ListQ[ColData],Length[ColData],1]}]
]
],
Null
]
]
},
{#2[[1]],ColI}
])&,YList]
)
]&,
data
],
Dimensions->{Length[YList],Length[data]},

FilterRules[{opts}, Options[Multipanel]],
XPlotRange->Map[Module[{ColData=#},(CalcPlotRange[If[ListQ[ColData],Flatten[Module[{YVar=#},#[All,{x,YVar}]&/@ColData]&/@YList],ColData[All,{x,#}]&/@YList]][[1]])]&,data],
YPlotRange->If[ListQ[data[[1]]],CalcPlotRange[Flatten@Module[{YVar=#},Map[#[All,{x,YVar}]&,data,{2}]]][[2]]&/@YList,CalcPlotRange[Module[{YVar=#},#[All,{x,YVar}]&/@data]][[2]]&/@YList],
Frame->True,
XFrameLabel->x,
YFrameLabel->YList,
YYShowTickLabels->Exterior,
YYShowFrameLabel->Exterior,
XXShowFrameLabel->If[OptionValue[XXFrameLabel]===Inherited,False,Exterior],

PanelLetter->None,
YPanelGaps->0.1,
XPanelGaps->0.1
],

FilterRules[{opts}, Options[Figure]],
CanvasSize->1*{6*Length[data],2*Length[YList]},
CanvasMargin->{{0.8,0.8},{0.6,0.6}},
Style->"DataAnalysis"
]
)
]


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


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