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

     htests.tl

     This file is part of the TREEFINDER software.

                                    COPYRIGHT (C) 1997-2008 BY GANGOLF JOBB
==============================================================================*)

SimulateSequences[
 hyp
][
 tree,model,nsites,prat,fnam,pkey,plen,gwel,nor
]:=<<

 hyp=LoadReport[hyp,AllSimulationData->True]|1,

 tree=hyp|'Phylogeny,
 model=hyp|'SubstitutionModel,
 nsites=hyp|'NSites,
 prat=hyp|'PartitionRates,

 If[

  IsSymbol[prat],

  GenerateSequences[tree,model,nsites],

  fnam=hyp|'FilterNames;
  pkey=hyp|'PartitionKeys;
  plen=hyp|'NSitesPartitionwise;
  gwel=hyp|'GroupwiseEdgeLengths;

  nor=
   If[IsSymbol[gwel],
    Dot[prat,plen/nsites],
    Table[Sum[i|1==g?i|2*j?0.,(=i,prat,j,plen=)]/Sum[i|1==g?j?0.,(=i,prat,j,plen=)],(=g,1:Max[prat]|1=)]
   ];

  IncludeFilters[
   Table[
    GenerateSequences[
     If[IsChain[j],
      ScaleEdgeLengths[gwel|(j|1),j|2/nor|(j|1)],
      ScaleEdgeLengths[tree,j/nor]
     ],
     If[IsList[model],model|k,model],
     i
    ],
    (=i,plen,j,prat,k,1:Length[pkey]=)
   ]@JoinSequences,
   fnam,pkey,plen
  ]

 ]

>>,
SetMask[SimulateSequences,"1"], Protect[SimulateSequences],

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

ParametricBootstrapTest[
 h0,
 h1,
 NReplicates->1000,
 Criterion->"Likelihood",
 ReturnPValue->False
][
 corr0,corr1,both,opt0,opt1,fil1:=~Void~,sim,rec0,rec1,lr:={},pval,msg
]:=<<

 both=JoinReports[h0,h1,AllHypothesisData->True],   (* loads reports and checks compatibility *)

 h0=both|1,
 h1=both|2,

 corr0=
  Criterion=="AIC"   ? (h0|'Likelihood)+0.5*(h0|'AIC )  ?
  Criterion=="AICc"  ? (h0|'Likelihood)+0.5*(h0|'AICc)  ?
  Criterion=="HQ"    ? (h0|'Likelihood)+0.5*(h0|'HQ)    ?
  Criterion=="BIC"   ? (h0|'Likelihood)+0.5*(h0|'BIC )  ?
  Criterion=="SimIC" ? (h0|'Likelihood)+0.5*(h0|'SimIC) ? 0.
 ,
 corr1=
  Criterion=="AIC"   ? (h1|'Likelihood)+0.5*(h1|'AIC )  ?
  Criterion=="AICc"  ? (h1|'Likelihood)+0.5*(h1|'AICc)  ?
  Criterion=="HQ"    ? (h1|'Likelihood)+0.5*(h1|'HQ)    ?
  Criterion=="BIC"   ? (h1|'Likelihood)+0.5*(h1|'BIC )  ?
  Criterion=="SimIC" ? (h1|'Likelihood)+0.5*(h1|'SimIC) ? 0.
 ,

 opt0={
  Tree->h0|'Phylogeny,
  SubstitutionModel->h0|'OSubstitutionModel,
  If[IsMember[h0,OPartitionRates->_],PartitionRates->h0|'OPartitionRates],
  If[IsMember[h0,OEdgeOptimizationOff->_],EdgeOptimizationOff->h0|'OEdgeOptimizationOff],
  AcceptFlatness->True,
  WithEdgeSupport->False
 },

 opt1={
  Tree->h1|'Phylogeny,
  SubstitutionModel->h1|'OSubstitutionModel,
  If[IsMember[h1,OPartitionRates->_],PartitionRates->h1|'OPartitionRates],
  If[IsMember[h1,OEdgeOptimizationOff->_],EdgeOptimizationOff->h1|'OEdgeOptimizationOff],
  AcceptFlatness->True,
  WithEdgeSupport->False
 },

 If[IsMember[h1,FilterNames->_],fil1={h1|'FilterNames,h1|'PartitionKeys,h1|'NSitesPartitionwise}],

 For[

  msg=" "<>i<>" ";Spell[msg];

  sim=SimulateSequences[h0];

  rec0=Prepend[opt0,sim]@ReconstructPhylogeny;

  If[!IsVoid[fil1],sim=Prepend[fil1,RemoveFilters[sim]]@IncludeFilters];

  rec1=Prepend[opt1,sim]@ReconstructPhylogeny;

  AppendR[lr,(rec1|1|'Likelihood-corr1)-(rec0|1|'Likelihood-corr0)];

  For[Spell["\b \b"],(=j,1:StringLength[msg]=)];

  Gc;

  ,(=i,1:NReplicates=)
 ],

 pval=CDF[SampleD[lr],(h1|'Likelihood-corr1)-(h0|'Likelihood-corr0)],

 If[ReturnPValue,Return[pval]],

 EditLines[{
   "",
   "PARAMETRIC BOOTSTRAP TEST",
   "=========================",
   "",
   "# Replicates: "<>NReplicates,
   "",
   "Criterion: "<>Criterion,
   "",
   "Simulated p-Value: "<>pval,
   ""
 }]

>>,
SetMask[ParametricBootstrapTest,"11."], Protect[ParametricBootstrapTest],
Create[ReturnPValue], Protect[ReturnPValue],

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

ComputeSimIC[
 model,
 NReplicates->1000
][
 mean,data,r0,r1
]:=<<

 model=
  LoadReport[model,AllHypothesisData->True]
   =|{(SimIC->_)->(),(NEffParameters->_)->()}
 ,

 Table[
  mean=0;
  For[
   data=SimulateSequences[{hyp}];
   r0={data,
       Tree->hyp|'Phylogeny,
       SubstitutionModel->hyp|'SubstitutionModel,
       If[IsMember[hyp,PartitionRates->_],PartitionRates->hyp|'PartitionRates],
       EdgeOptimizationOff->True,
       WithEdgeSupport->False
      }@ReconstructPhylogeny;
   r1={data,
       Tree->hyp|'Phylogeny,
       SubstitutionModel->hyp|'OSubstitutionModel,
       If[IsMember[hyp,OPartitionRates->_],PartitionRates->hyp|'OPartitionRates],
       If[IsMember[hyp,OEdgeOptimizationOff->_],EdgeOptimizationOff->hyp|'OEdgeOptimizationOff],
       AcceptFlatness->True,
       WithEdgeSupport->False
      }@ReconstructPhylogeny;
   mean+=(r1|1|'Likelihood)-(r0|1|'Likelihood);
   Gc;
   ,(=i,1:NReplicates=)
  ];
  mean/=NReplicates;
  Join[{SimIC->-2*(hyp|'Likelihood-mean),NEffParameters->mean},hyp]
  ,(=hyp,model=)
 ]

>>,
SetMask[ComputeSimIC,"1."], Protect[ComputeSimIC],

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

ComputeMLDistances[
 data,
 hyp
][
 list:={}
]:=<<

 data=LoadSequences[data],

 hyp=LoadReport[hyp,AllHypothesisData->True,Sequences->'data]|1,

 {
  If[
   IsMember[hyp,FilterNames->_],
   SelectFilters[data,hyp|'FilterNames],
   data
  ],
  Tree->hyp|'Phylogeny, (* prevents reestimation *)
  SubstitutionModel->hyp|'SubstitutionModel,
  If[IsMember[hyp,PartitionRates->_],PartitionRates->hyp|'PartitionRates],
  If[IsMember[hyp,GroupwiseEdgeLengths->_],GroupwiseEdgeLengths->hyp|'GroupwiseEdgeLengths],
  EdgeOptimizationOff->True,
  WithEdgeSupport->False,
  AppendPairwiseDistances->'list
 }@ReconstructPhylogeny,

 Drop,

 TakeFirst[list]

>>,
SetMask[ComputeMLDistances,"11"], Protect[ComputeMLDistances],

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

ComputeSitewiseRates[
 data,
 hyp
][
 list:={}
]:=<<

 data=LoadSequences[data],

 hyp=LoadReport[hyp,AllHypothesisData->True,Sequences->'data]|1,

 {
  If[
   IsMember[hyp,FilterNames->_],
   SelectFilters[data,hyp|'FilterNames],
   data
  ],
  Tree->hyp|'Phylogeny,
  SubstitutionModel->hyp|'SubstitutionModel,
  If[IsMember[hyp,PartitionRates->_],PartitionRates->hyp|'PartitionRates],
  If[IsMember[hyp,GroupwiseEdgeLengths->_],GroupwiseEdgeLengths->hyp|'GroupwiseEdgeLengths],
  EdgeOptimizationOff->True,
  WithEdgeSupport->False,
  AppendSitewiseRates->'list
 }@ReconstructPhylogeny,

 Drop,

 TakeFirst[list]

>>,
SetMask[ComputeSitewiseRates,"11"], Protect[ComputeSitewiseRates],

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

ComputeSitewiseLikelihoods[
 data,
 model,
 Criterion->"Likelihood"
][
 list:={}
]:=<<

 data=LoadSequences[data],

 model=LoadReport[model,AllHypothesisData->True,Sequences->'data],

 For[
  {
   If[
    IsMember[hyp,FilterNames->_],
    SelectFilters[data,hyp|'FilterNames],
    data
   ],
   Tree->hyp|'Phylogeny,
   SubstitutionModel->hyp|'SubstitutionModel,
   If[IsMember[hyp,PartitionRates->_],PartitionRates->hyp|'PartitionRates],
   If[IsMember[hyp,GroupwiseEdgeLengths->_],GroupwiseEdgeLengths->hyp|'GroupwiseEdgeLengths],
   EdgeOptimizationOff->True,
   WithEdgeSupport->False,
   Verbose->False,
   AppendSitewiseLikelihoods->'list
  }@ReconstructPhylogeny,
  (=hyp,model=)
 ],

 TakeAll[list],

 Criterion!="Likelihood",
 <<
  -Table[((hyp|'Likelihood)+0.5*(hyp|Interpret[Criterion]))/(hyp|'NSites),(=hyp,model=)],
  'AddArrays,Map
 >>,If1

>>,
SetMask[ComputeSitewiseLikelihoods,"11."], Protect[ComputeSitewiseLikelihoods],
Create[Criterion], Protect[Criterion],

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

TestTopologies[
 data,
 trees,
 SubstitutionModel->~Void~,
 PartitionRates->~Void~,
 EdgeOptimizationOff->~Void~,
 NReplicates->50000
][
 list:={}
]:=<<

 data=LoadSequences[data],

 trees=LoadTreeList[trees],

 trees=
  Table[
   {
    data,
    Tree->hyp,
    If[!IsVoid[SubstitutionModel],SubstitutionModel->SubstitutionModel],
    If[!IsVoid[PartitionRates],PartitionRates->PartitionRates],
    If[!IsVoid[EdgeOptimizationOff],EdgeOptimizationOff->EdgeOptimizationOff],
    WithEdgeSupport->False,
    AppendSitewiseLikelihoods->'list
   }@ReconstructPhylogeny,
   (=hyp,trees=)
  ]
  @Join
 ,

 PairedSitesTests[TakeAll[list],NReplicates->NReplicates]#trees#Join

>>,
SetMask[TestTopologies,"11."], Protect[TestTopologies],

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

TestHypotheses[
 data,
 model,
 Criterion->"AICc",
 NReplicates->50000
]:=<<

 data=LoadSequences[data],

 model=
  LoadReport[model,AllHypothesisData->True,Sequences->'data]
   =|{(PairedSitesSupport->_)->(),(PairedSitesCriterion->_)->()}
 ,

 PairedSitesTests[
  ComputeSitewiseLikelihoods[data,model,Criterion->Criterion],
  NReplicates->NReplicates
 ]
 #Table[{PairedSitesCriterion->Criterion},(=i,model=)]
 #model
 #Join

>>,
SetMask[TestHypotheses,"11."], Protect[TestHypotheses],

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

ProposeModel[
 data,
 Criterion->"AICc",
 ModelType->"DNA",
 Tree->~Void~,
 PlusF->True,
 PlusG->True,
 PlusI->True,
 PlusGI->True,
 WithOptimizedFrequencies->False,
 NRateCategories->5,
 Partitionwise->False,
 Candidates->~Void~,
 BasicCandidates->~Void~,
 ShowCandidates->False,
 ShowParameters->False,
 ShowReport->False,
 ReturnReport->False
][
 result,msg
]:=<<

 data=LoadSequences[data],

 If[IsVoid[Candidates],
  Candidates=
   ModelType=="DNA"     ? '{HKY[{3,1,1,1,1,3},Empirical],HKY[Optimum,Empirical],TN[Optimum,Empirical],J1[Optimum,Empirical],J2[Optimum,Empirical],J3[Optimum,Empirical],TVM[Optimum,Empirical],GTR[Optimum,Empirical]} ?
   ModelType=="DNA*"    ? '{HKY[{3,1,1,1,1,3},Empirical],HKY[Optimum,Empirical],TN[Optimum,Empirical],GTR[Optimum,Empirical]} ?
   ModelType=="Protein" ? '{Dayhoff[,],BLOSUM[,],JTT[,],LG[,],PMB[,],VT[,],WAG[,],mtREV[,],mtMam[,],mtArt[,],cpREV[,],rtREV[,],betHIV[,],witHIV[,],MIX[Dayhoff,BLOSUM,JTT,LG,PMB,VT,WAG,mtREV,mtMam,mtArt,cpREV,rtREV,betHIV,witHIV][Optimum,]} ?
   ModelType=="NuProt"  ? '{Dayhoff[,],BLOSUM[,],JTT[,],LG[,],PMB[,],VT[,],WAG[,],MIX[Dayhoff,BLOSUM,JTT,LG,PMB,VT,WAG][Optimum,]} ?
   ModelType=="MtProt"  ? '{mtREV[,],mtMam[,],mtArt[,],MIX[mtREV,mtMam,mtArt][Optimum,]} ?
   ModelType=="CpProt"  ? '{cpREV[,]} ?
   ModelType=="RtProt"  ? '{rtREV[,],betHIV[,],witHIV[,],MIX[rtREV,betHIV,witHIV][Optimum,]} ?
   ModelType=="RNA"     ? '{bactRNA[,],eukRNA[,],euk23RNA[,],mitoRNA[,],MIX[bactRNA,eukRNA,euk23RNA,mitoRNA][Optimum,]} ?
   ModelType=="DG"      ? '{DG[Optimum,Empirical]} ?
   ModelType=="DNA3"    ? '{GTR3[Optimum,Empirical]} ?
   ModelType=="DNA2"    ? '{GTR2[{},Empirical]} ?
   ModelType=="UserDef" ? ( IsList[BasicCandidates]&&BasicCandidates#IsSubstitutionModel@And ? Members[BasicCandidates^|{_s:__->_s}] ? Throw["ProposeModel: Bad user-defined candidates."] ) ?
   Throw["ProposeModel: Undefined data type."];
  If[ PlusF                    , Candidates=Candidates^|{ _m[_r,]->(_m[_r,],_m[_r,Empirical]) } ];
  If[ WithOptimizedFrequencies , Candidates=Candidates^|{ _m[_r,Empirical]->_m[_r,Optimum]    } ];
  Candidates=
   Join[         Candidates                                                         ,
    If[ PlusG  , Candidates^|`{ _m[_r,_f]->_m[_r,_f]:G [Optimum]:NRateCategories } ],
    If[ PlusI  , Candidates^|`{ _m[_r,_f]->_m[_r,_f]:I [Optimum]                 } ],
    If[ PlusGI , Candidates^|`{ _m[_r,_f]->_m[_r,_f]:GI[Optimum]:NRateCategories } ]
   ]
 ],

 If[ShowCandidates,EditList[Candidates,Description->"\n(* set of candidate models *)\n\n"]],

 If[IsVoid[Tree],
  Tree=
   ReconstructPhylogeny[
    data,
    SubstitutionModel->Candidates|1-|{ _m[_r,Optimum]->_m[_r,Empirical] },
    SearchDepth->1,
    AcceptFlatness->True,
    WithEdgeSupport->False
   ]|1|'Phylogeny
 ],

 If[Partitionwise,

  If[!HasFilters[data],Throw["ProposeModel: Data is expected to be partitioned."]];

  result=
   Table[
    msg=" #"<>k<>" ";Spell[msg];
    ProposeModel[
     SelectPartitions[data,{k}],
     Criterion->Criterion,
     Tree->Tree,
     Partitionwise->False,
     Candidates->Candidates,
     ReturnReport->True
    ]#<<{(Meaning->_str)=>(Meaning->"Partition #"<>k<>" / "<>_str)},ReplaceContents>>
    ;For[Spell["\b \b"],(=j,1:StringLength[msg]=)]
    ,(=k,GetPartitionKeys[data]=)
   ];

  If[ReturnReport,Return[result@Join]];

  If[ShowReport,PlotTree[result@Join]];

  If[ShowParameters,EditList[result#<<1,Get,'SubstitutionModel,Get>>,Description->"\n(* estimated parameters of partitionwise best-fit models *)\n\n"]];

  Return[result#<<1,Get,'OSubstitutionModel,Get>>]

 ],

 result=
  SortHypotheses[
   Table[
    msg=" "<>i<>"/"<>Length[Candidates]<>" ";Spell[msg];
    ReconstructPhylogeny[
     data,
     Tree->Tree,
     SubstitutionModel->m,
     AcceptFlatness->True,
     WithEdgeSupport->False
    ]|1
    ;For[Spell["\b \b"],(=j,1:StringLength[msg]=)]
    ,(=m,Candidates,i,1:Length[Candidates]=)
   ],
   If[Criterion=="Likelihood",Criterion,"-"<>Criterion]
  ],

 result=Table[Prepend[h,Meaning->Criterion<>" Rank "<>i],(=h,result,i,1:Length[result]=)];

 If[ReturnReport,Return[result]],

 If[ShowReport,PlotTree[result]],

 If[ShowParameters,Edit[result|1|'SubstitutionModel,Description->"\n(* estimated parameters of best-fit model *)\n\n"]],

 result|1|'OSubstitutionModel

>>,
SetMask[ProposeModel,"1."], Protect[ProposeModel],
Create[ModelType], Protect[ModelType],
Create[PlusF], Protect[PlusF],
Create[PlusG], Protect[PlusG],
Create[PlusI], Protect[PlusI],
Create[PlusGI], Protect[PlusGI],
Create[WithOptimizedFrequencies], Protect[WithOptimizedFrequencies],
Create[NRateCategories], Protect[NRateCategories],
Create[Partitionwise], Protect[Partitionwise],
Create[Candidates], Protect[Candidates],
Create[BasicCandidates], Protect[BasicCandidates],
Create[ShowCandidates], Protect[ShowCandidates],
Create[ShowParameters], Protect[ShowParameters],
Create[ShowReport], Protect[ShowReport],
Create[ReturnReport], Protect[ReturnReport],

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

()
