(*^ ::[ 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[] (* ^*)