(*^

::[	frontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.1";
	macintoshStandardFontEncoding; 
	currentKernel; 
	fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, "Times"; ;
	fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, "Times"; ;
	fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, "Times"; ;
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, "Times"; ;
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, "Times"; ;
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, "Times"; ;
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; ;
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5,  12, "Courier"; ;
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; ;
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5,  12, "Courier"; ;
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; ;
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5,  12, "Courier"; ;
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287,  12, "Courier"; ;
	fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  10, "Geneva"; ;
	fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = leftheader, inactive, L2,  12, "Times"; ;
	fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7,  12, "Times"; ;
	fontset = leftfooter, inactive, L2,  12, "Times"; ;
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; ;
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; ;
]
:[font = input; initialization; preserveAspect; startGroup; ]
*)
BeginPackage["DDsim`"]
(*
:[font = output; output; inactive; preserveAspect; endGroup; ]
"DDsim`"
;[o]
DDsim`
:[font = input; inactive; preserveAspect; ]
Copyright  A.Grafen and M. Ridley   1997
:[font = input; inactive; preserveAspect; ]
(* Need to add a parameter checking level to the trial *)
:[font = input; initialization; preserveAspect; startGroup; ]
*)
trial::usage="Produces simulated discrete datasets following\n "<>
"Grafen and Ridley (1997. J. theor. Biol. 184: 7-14.).\nParameters are\n"<>
" 1. name of output file\n"<>
" 2. a list with ten elements:\n"<>
"    (i) number of simulated datasets required\n"<>
"   (ii) initial states for C and D, as {c0,d0}\n"<>
"  (iii) the value of rho used to transform node heights\n"<>
"   (iv) the rate of evolution\n"<>
"    (v) the phylogeny\n"<>
"   (vi) rule for determining node heights\n"<>
"  (vii) transition-rate matrix for C\n"<>
" (viii) transition-rate matrix for D\n"<>
"   (ix) rule for determining A from B and C\n"<>
"    (x) rule for determining B from D\n\n"<>
"See accompanying file 'Example use' for worked examples.";
(*
:[font = output; output; inactive; initialization; preserveAspect; fontLeading = 0; ]
"Produces simulated discrete datasets following Grafen and \nRidley \
  (1997. A new model for discrete character evolution.\nJ. theor. Biol. \
  184, 7-14). Parameters are\n 1. name of output file\n 2. a list with nine \
  elements:\n  (i) number of simulated datasets required\n (ii) initial \
  states for C and D, as {c0,d0}\n (iii) the value of rho used to transform \
  node heights\n (iv) the rate of evolution\n (v) the phylogeny\n (v) rule \
  for determining node heights\n (vii) transition-rate matrix for C\n \
  (viii) transition-rate matrix for D\n (ix) rule for determining A from \
  C\n (ix) rule for determining B from A and D\n\nSee accompanying file \
  'Example use' for worked examples."
;[o]
Produces simulated discrete datasets following Grafen and 
Ridley (1997. A new model for discrete character evolution.
J. theor. Biol. 184, 7-14). Parameters are
 1. name of output file
 2. a list with nine elements:
  (i) number of simulated datasets required
 (ii) initial states for C and D, as {c0,d0}
 (iii) the value of rho used to transform node heights
 (iv) the rate of evolution
 (v) the phylogeny
 (v) rule for determining node heights
 (vii) transition-rate matrix for C
 (viii) transition-rate matrix for D
 (ix) rule for determining A from C
 (ix) rule for determining B from A and D

See accompanying file 'Example use' for worked examples.
:[font = input; initialization; preserveAspect; endGroup; ]
*)
symtree::usage="Generates symmetrical trees with #1 twigs per"<>
" higher node, and #2 levels below the root, leading to a"<>
" tree with #1^#2 species";
(*
:[font = input; initialization; preserveAspect; startGroup; ]
*)
Begin["`Private`"]
(*
:[font = output; output; inactive; preserveAspect; endGroup; ]
"DDsim`Private`"
;[o]
DDsim`Private`
:[font = input; dontPreserveAspect; ]
(*This is based on pfa 1.1, the apparently most modern
  of the pfa files.  It differs in having maptreeb from
  pfa 1.0 (and that was an essential correction to
  earlier versions of pfa, which had errors in the
  use of rho)
  *)
:[font = input; preserveAspect; ]
(*I have hard-wired treat as the form_, for simplicity of
  user interface, and omitted the first argument of treat
  from definition and use.  (The first argument wasn't
  used in production runs anyway - it may once have been
  used for debugging.) *)
:[font = input; preserveAspect; ]
(* I have un-hardwired arule and brule, to allow them to
   be specified by the user. *)
:[font = input; initialization; preserveAspect; ]
*)
trial[outfile_String, parmlist_List] := 
    If[check[parmlist],
      (simsout = OpenWrite[outfile]; Write[simsout,parmlist];
       Print[Timing[Apply[dosims, convert[parmlist]]]];
       Close[simsout]),
      Print["Simulation not performed because parameters unsuitable"]];

check[l_List]:=True

convert[l_List]:=
  MapAt[
    ReplaceAll[#,{Global`tip->DDsim`Private`tip}]&,
    l,
    5]   
    
 
dosims[nsims_, {c_, d_}, rho_, rate_, tree_, htf_, ctrans_, dtrans_, 
   arule_, brule_] := (ap = allpairs[tree]; holderc =. ; 
    newmearray[holderc, MatrixExp[ctrans], Last[ap][[1]], rho, rate, htf, 
     ap]; holderd =. ; newmearray[holderd, MatrixExp[dtrans], Last[ap][[1]], 
     rho, rate, htf, ap]; (realisma[preparecomrefma[tree], holderc, 
        holderd, c, d, arule, brule] & ) /@ Range[1, nsims])
 
(* ctrans = {{-0.01, 0.01, 0., 0., 0., 0.}, {0.2, -0.82, 0.6, 0.02, 0., 0.}, 
   {0., 0.6, -0.6199999999999999999, 0.02, 0., 0.}, 
   {0., 0., 0.02, -0.6199999999999999999, 0.6, 0.}, 
   {0., 0., 0.02, 0.6, -0.82, 0.2}, {0., 0., 0., 0., 0.01, -0.01}} *)
 
(* dtrans = {{-0.01, 0.01, 0., 0., 0., 0.}, {0.2, -0.82, 0.6, 0.02, 0., 0.}, 
   {0., 0.6, -0.6199999999999999999, 0.02, 0., 0.}, 
   {0., 0., 0.02, -0.6199999999999999999, 0.6, 0.}, 
   {0., 0., 0.02, 0.6, -0.82, 0.2}, {0., 0., 0., 0., 0.01, -0.01}} *)
 
allpairs[tree_List] := 
 Union[Flatten[Select[Flatten[maptree[addtopairs[#1] & , Null & , Null & , 
       addtipcount[tree, Identity], dummy]], Head[#1] === pt & ], 3, pt]]

maptree[f_, g_, h_, tree_List, carry_] := 
  maptreeb[f, g, h, tree, carry, g[tree, carry]]
 
maptree[f_, g_, h_, tree_tip, carry_] := h[tree, carry]
 
maptreeb[f_, g_, h_, tree_List, carry_, newcarry_] := 
  Prepend[(maptree[f, g, h, #1, newcarry] & ) /@ Rest[tree], f[tree, carry]]
 
addtopairs[tree_List] := Apply[pt, pairs[First /@ Rest[tree]]]
 
pairs[l_List] := Union[Flatten[fixtop /@ subsetsbetw[l, 2, Length[l]], 1]]
 
fixtop[l_List] := 
  ({Apply[Plus, l], Apply[Plus, #1]} & ) /@ subsetsbetw[l, 1, Length[l] - 1]
 
subsetsbetw[l_List, i_Integer, j_Integer] := 
  Select[allsubsets[l, j], Length[#1] >= i & ]
 
allsubsets[{a_}] := {{}, {a}}
 
allsubsets[l_List, k_Integer] := 
  Apply[Join, Table[allsubsetsfs[l, jj], {jj, 0, k}]]
 
allsubsetsfs[{a___}, 0] := {{}}
 
allsubsetsfs[{a_}, 1] := {{a}}
 
allsubsetsfs[l_List, j_Integer] := 
  If[j == Length[l], {l}, Join[(Prepend[#1, l[[1]]] & ) /@ 
     allsubsetsfs[Rest[l], j - 1], allsubsetsfs[Rest[l], j]]]
 
addtipcount[l_List, op_] := (countme[#1, op] & ) //@ l
 
countme[t_tip, op_] := Prepend[t, op[1]]
 
countme[l_List, op_] := Prepend[Rest[l], Apply[Plus, First /@ Rest[l]]]
 
countme[i_Integer, op_] := i
 
newmearray[name_, trans_, n_, rho_, rate_, htf_, pairs_] := 
  Block[{i, j}, name =. ; name = 
     Table[Table[If[MemberQ[pairs, {j, i}], 
        parsums /@ 
         MatrixPower[trans, rate*
           N[(htf[j]/htf[n])^rho - (htf[i]/htf[n])^rho]], "blank cell"], 
       {i, 1, j}], {j, 1, n}]; Null]
 
parsums[l_] := Block[{temp = 0}, ((temp = temp + #1) & ) /@ l]
 
realisma[tree_, m_List, n_List, c_Integer, d_Integer, arule_,
         brule_] := 
  treat[abrealis[cdrealisma[tree, m, n, c, d],arule,brule]]
 
abrealis[tree_,arule_,brule_] :=
        hittips[tip[#1[[1]], speciesdata[#1[[2]],arule,brule]] & , tree]
 
hittips[op_, tree_List] := 
  Prepend[(hittips[op, #1] & ) /@ Rest[tree], First[tree]]
 
hittips[op_, tree_tip] := op[tree]
 
(* I've swapped ci and di round in the second line here, in the
   hope of compatibility with the MS *) 
speciesdata[{ci_Integer, di_Integer},arule_,brule_] :=
    ({arule[#1, ci], #1} & )[brule[di]]
 
(* arule[bi_Integer, di_Integer] := {1, 1, 2, 1, 2, 2}[[di]] *)
 
(* brule[ci_Integer] := {1, 1, 2, 1, 2, 2}[[ci]] *)
 
cdrealisma[tree_, m_List, n_List, c_Integer, d_Integer] := 
  delaymaptree[#1[[1]] & , {nextvaluema[#2[[1]], m, #2[[3]], #1[[1]]], 
     nextvaluema[#2[[2]], n, #2[[3]], #1[[1]]], #1[[1]]} & , 
   Append[#1, {nextvaluema[#2[[1]], m, #2[[3]], 1], 
      nextvaluema[#2[[2]], n, #2[[3]], 1]}] & , tree, {c, d, tree[[1]]}]
 
delaymaptree[f_, g_, h_, tree_List, carry_] := 
  Prepend[(maptree[f, g, h, #1, carry] & ) /@ Rest[tree], First[tree]]
 
nextvaluema[ini_Integer, m_List, i_Integer, j_Integer] := 
  vindex[m[[i,j,ini]] - Random[]]
 
vindex[l_] := Position[l, _Real?NonNegative][[1,1]]
 
preparecomrefma[tree_] := addnodeheights[refine[numbertips[tree]], Identity]
 
addnodeheights[l_List, op_] := 
  maptree[#1[[1]] & , Null & , Drop[#1, 1] & , addtipcount[l, op], dummy]
 
refine[l_List] := Prepend[refinenode[Rest[l], refine], First[l]]
 
refine[l_tip] := l
 
refinenode[l_List, op_] := putparts1[l, ransplit[Length[l]], op]
 
refinenode[l_tip, op_] := l
 
putparts1[l_List, tree_List, op_] := (putparts2[l, #1, op] & ) /@ tree
 
putparts1[l_List, i_Integer, op_] := op[l[[i]]]
 
putparts2[l_List, tree_List, op_] := 
  Prepend[(putparts2[l, #1, op] & ) /@ tree, Null]
 
putparts2[l_List, i_Integer, op_] := op[l[[i]]]
 
ransplit[i_Integer] := split[ranreorderrange[i]]
 
split[(l_List)?(Length[#1] >= 3 & )] := 
  split /@ ({Take[l, #1], Drop[l, #1]} & )[ransplitn[Length[l]]]
 
split[(l_List)?(Length[#1] == 2 & )] := l
 
split[(l_List)?(Length[#1] == 1 & )] := l[[1]]
 
ransplitn[i_Integer] := 1 + Floor[(i - 1)*Random[]]
 
ranreorderrange[i_] := Last /@ Sort[({Random[], #1} & ) /@ Range[i]]
 
numbertips[tree_List] := 
  Block[{ii = 1}, (If[Head[#1] === tip, tip[ii++], #1] & ) //@ tree]
  
treat[x_] := 
  (Write[simsout, gettips[x]];)

gettips[l_List] := 
  (#1[[2]] & ) /@ Sort[Select[Flatten[l], Head[#1] === tip & ]]

symtree[twigs_, 1] := Array[If[#1 == 1, Null, Global`tip[]] & , {twigs + 1}]
 
symtree[twigs_, depth_] := 
  Array[If[#1 == 1, Null, symtree[twigs, depth - 1]] & , {twigs + 1}]
 

(*
:[font = input; initialization; preserveAspect; startGroup; ]
*)
End[]
(*
:[font = output; output; inactive; preserveAspect; endGroup; ]
"DDsim`Private`"
;[o]
DDsim`Private`
:[font = input; initialization; preserveAspect; ]
*)
EndPackage[]
(*
^*)