
(*==============================================================================

     samplestat.tl

     This file is part of the TREEFINDER software.


                                    COPYRIGHT (C) 1997-2008 BY GANGOLF JOBB

==============================================================================*)

ThreadQuantile[data,p][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadQuantile: Bad data."]],
 If[p<0||p>1,Throw["ThreadQuantile: Second argument is expected to be a probability."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadQuantile: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#<<'d,Def,Quantile[SampleD[d],p],'d,Forget>>]
>>,
SetMask[ThreadQuantile,"11"], Protect[ThreadQuantile],

(*============================================================================*)

ThreadMean[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadMean: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadMean: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#Mean]
>>,
SetMask[ThreadMean,"1"], Protect[ThreadMean],

(*============================================================================*)

ThreadMeanInterval[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadMeanInterval: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadMeanInterval: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#<<'d,Def,Mean[d][Quantile[SampleD[d],0.025],Quantile[SampleD[d],0.975]],'d,Forget>>]
>>,
SetMask[ThreadMeanInterval,"1"], Protect[ThreadMeanInterval],

(*============================================================================*)

ThreadMedian[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadMedian: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadMedian: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#Median]
>>,
SetMask[ThreadMedian,"1"], Protect[ThreadMedian],

(*============================================================================*)

ThreadMaximum[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadMaximum: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadMaximum: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#Max]
>>,
SetMask[ThreadMaximum,"1"], Protect[ThreadMaximum],

(*============================================================================*)

ThreadMinimum[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadMinimum: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadMinimum: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#Min]
>>,
SetMask[ThreadMinimum,"1"], Protect[ThreadMinimum],

(*============================================================================*)

ThreadRange[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadRange: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadRange: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#Range]
>>,
SetMask[ThreadRange,"1"], Protect[ThreadRange],

(*============================================================================*)

ThreadVariance[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadVariance: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadVariance: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#Variance]
>>,
SetMask[ThreadVariance,"1"], Protect[ThreadVariance],

(*============================================================================*)

ThreadStandardDeviation[data][first]:=<<
 If[!IsList[data]||Length[data]<2,Throw["ThreadStandardDeviation: Bad data."]],
 first=ZeroNumbers[data|1],For[If[ZeroNumbers[i]!=first,Throw["ThreadStandardDeviation: Incompatible sample."]],(=i,data=)],
 PutNumbers[first,Transpose[data#GetNumbers]#StandardDeviation]
>>,
SetMask[ThreadStandardDeviation,"1"], Protect[ThreadStandardDeviation],

(*============================================================================*)

CalibrateTreeRaw:=Lookup[CalibrateTree], Protect[CalibrateTreeRaw],

CalibrateTree[
 data,
 calib,
 Method->~Void~,
 MinimumEdgeLength->~Void~,
 NTimeReplicates->0,
 SampleFile->~Void~
][
 fun,sf
]:=<<
 data=AsTreeList[data];
 calib=LoadCalibrationData[calib];
 fun=<<
  calib,
  `If[!IsVoid[Method],Method->Method],
  `If[!IsVoid[MinimumEdgeLength],MinimumEdgeLength->MinimumEdgeLength],
  CalibrateTreeRaw
 >>;
 data=
  If[NTimeReplicates>0,
   WrapR[fun,1,'ReplaceContents[{_p:_a:_b=>_p:RandomReal[_a,_b]}]];
   Table[data#fun,(=i,1:NTimeReplicates=)]@Join,
   TakeAll[data]#fun
  ];
 If[!IsVoid[SampleFile],
  sf=Open[SampleFile,"w"];Spell["{\n",sf];
  data#<<" ",sf,Spell,sf,Write,",\n",sf,Spell>>;
  Spell[" ()\n}\n",sf];Close[sf]
 ];
 If[Length[data]>1,data={ThreadMeanInterval[data]}];
 data|1
>>,
SetMask[CalibrateTree,"11."], Protect[CalibrateTree],
Create[NTimeReplicates], Protect[NTimeReplicates],

(*============================================================================*)

()
