︠addb8005-1594-4574-9721-90f6c6cdc969i︠ %html
This is code accompanying the technical report 'The Transmission Process' by Raazesh Sainudiin and David Welch. The code is for clarity and simplicity of implementation (it is not efficient for large scale simulations).
︡bcca25c8-9302-4f45-ae5b-69ff475eb56e︡{"html":"\nThis is code accompanying the technical report 'The Transmission Process' by Raazesh Sainudiin and David Welch. The code is for clarity and simplicity of implementation (it is not efficient for large scale simulations).\n
\n"}︡ ︠c9e0dc04-2964-47ac-ae31-caa24e41ad71i︠ ︡355398ba-813a-4eed-9309-e4a6b8c29326︡ ︠3eb3fe29-0afd-4fd5-a418-e3127eaec3f3s︠ def CountsDict(X): '''convert a list X into a Dictionary of counts or frequencies''' CD = {} for x in X: CD[x] = (CD[x] + 1) if (x in CD) else 1 return CD def markAsInfected(C,v,m): '''mark node v as infected with marker m on each of the incoming edges of v in SICN C''' for e in C.incoming_edge_iterator([v]): C.set_edge_label(e[0],e[1],m) def susceptibleOutEdges(C,vs): '''return the the susceptible outedges of node v in vs in SICN C''' SOE = [e for e in C.outgoing_edge_iterator(vs) if e[2]==None] return SOE def growTransmissionTree(Ttree, pDict, z, infector, infectee): '''grow the transmission tree Ttree and update pathsDict pDict by adding the z-th infection event with infector -> infectee ''' LBT = LabelledBinaryTree newSubTree = LBT([LBT([None,None], label=infector), LBT([None, None], label=infectee)], label=z).clone() path2Infector = pDict[infector] if z==1: Ttree = newSubTree else: Ttree[tuple(path2Infector)] =newSubTree #print ascii_art(Ttree) pDict[infector]=path2Infector+[0] pDict[infectee]=path2Infector+[1] pDict[z]=path2Infector return Ttree def forgetLeafLabels(T): '''return the transmission tree T with all leaf labels set to 0''' leafLabelSet=set(T.leaf_labels()) leafUnlabelledT=T.map_labels(lambda z:0 if (z in leafLabelSet) else z) return leafUnlabelledT def forgetAllLabels(T): '''return the transmission tree T with all node labels removed''' return T.shape() def justTree(T): '''return the transmission tree T as nonplanar unranked unlabelled tree''' return Graph(T.shape().to_undirected_graph(),immutable=True) def transmissionProcessTC(C,initialI): '''return transmission tree outcome of the DTDS transmission MC on SICN C with initial infection at node initialI''' #initialisation of SICN z=0 # infection event count markAsInfected(C,initialI,'infected') infectedIs = [initialI] popSize=C.order() # initialisation of Transmission Tree pathsDict={} # dictionary of nodes -> paths from root in tree LBT = LabelledBinaryTree # individuals in tree are labelled by "i"+str(integer_label) T = LBT([None,None],label="i"+str(initialI)).clone() pathsDict["i"+str(initialI)]=[] while (len(infectedIs) < popSize): z=z+1 # increment infection event count currentSOE = susceptibleOutEdges(C,infectedIs) numberInfected=len(currentSOE) nextEdge = currentSOE[randrange(0,numberInfected)] C.set_edge_label(nextEdge[0],nextEdge[1],z) infectedIs.append(nextEdge[1]) markAsInfected(C,nextEdge[1],'inf') T=growTransmissionTree(T, pathsDict, z, "i"+str(nextEdge[0]),"i"+str(nextEdge[1])) #print "step z = ",z; print ascii_art(T); print "--------------------" return T.as_ordered_tree(with_leaves=False) ︡df385f0b-32e9-4477-9ef9-3d0473a49788︡{"done":true} ︠61f20284-ea37-4ddf-aa0c-914903624397i︠ %md # Three example host contact networks ︡c7a32ce2-f758-4cc7-ab0d-87d122fd26d0︡{"done":true,"md":"# Three example host contact networks"} ︠1328509d-e2f8-4890-ba8e-035c95242c78︠ ︡08006af5-7dcd-4f61-97d2-ddd38b605a7c︡ ︠6626bde8-e783-49b1-ab18-de1ccc1bbca7i︠ %htmlLet the SICN be the complete graph. We get the distribution of transmission trees at various resolutions next.
︡8c2b1377-5069-4252-a12d-f37a636bc4ed︡{"html":"\n Let the SICN be the complete graph. We get the distribution of transmission trees at various resolutions next.\n
\n"}︡ ︠c8050ada-0aa4-4cf3-8721-a207dd1c45ffs︠ graphs.CompleteGraph(4).to_directed().show() ︡73c63af1-d4aa-4efd-9761-c4b76202cb25︡{"file":{"filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute5-us/12878/tmp_XdZbf1.svg","show":true,"text":null,"uuid":"5f073793-f0eb-4b3d-87a4-e424867b553c"},"once":false}︡{"done":true}︡ ︠fdcc0361-a170-4edf-a2b6-49421e597e14︠ # demo transmissionProcessTC(graphs.CompleteGraph(4).to_directed(),0) ︡088fd35a-32e2-4d56-b17f-4f4aa047080f︡{"stdout":"step z = "}︡{"stdout":" 1\n 1_\n / \\\ni0 i3"}︡{"stdout":"\n--------------------\nstep z = 2\n __1__\n / \\\n 2_ i3\n / \\ \ni0 i2 \n--------------------\nstep z = "}︡{"stdout":" "}︡{"stdout":"3"}︡{"stdout":"\n"}︡{"stdout":" ____1____\n / \\\n _2__ i3\n / \\ \ni0 3_ \n / \\ \n i2 i1 "}︡{"stdout":"\n--------------------\n1[2[i0[], 3[i2[], i1[]]], i3[]]"}︡{"stdout":"\n"}︡{"done":true}︡ ︠31eb74e5-db23-4729-aac5-937d5c5e5c96︠ ts=[transmissionProcessTC(graphs.CompleteGraph(4).to_directed(),0) for _ in range(100000)] d=CountsDict(ts) print len(d) d ︡cffa9b82-ebaf-4224-ba7d-609251b50d93︡{"stdout":"36\n"}︡{"stdout":"{1[3[i0[], i2[]], 2[i1[], i3[]]]: 2824, 1[i0[], 2[i1[], 3[i3[], i2[]]]]: 2752, 1[i0[], 2[3[i1[], i3[]], i2[]]]: 2799, 1[2[i0[], i2[]], 3[i1[], i3[]]]: 2793, 1[2[i0[], 3[i1[], i2[]]], i3[]]: 2650, 1[2[i0[], 3[i1[], i3[]]], i2[]]: 2717, 1[3[i0[], i3[]], 2[i2[], i1[]]]: 2726, 1[i0[], 2[3[i3[], i1[]], i2[]]]: 2780, 1[2[i0[], i3[]], 3[i2[], i1[]]]: 2771, 1[i0[], 2[3[i3[], i2[]], i1[]]]: 2762, 1[2[3[i0[], i1[]], i3[]], i2[]]: 2770, 1[2[i0[], 3[i2[], i3[]]], i1[]]: 2751, 1[i0[], 2[i2[], 3[i1[], i3[]]]]: 2850, 1[i0[], 2[i3[], 3[i1[], i2[]]]]: 2719, 1[2[i0[], 3[i2[], i1[]]], i3[]]: 2835, 1[3[i0[], i2[]], 2[i3[], i1[]]]: 2793, 1[3[i0[], i1[]], 2[i3[], i2[]]]: 2823, 1[2[i0[], i1[]], 3[i2[], i3[]]]: 2815, 1[2[i0[], i1[]], 3[i3[], i2[]]]: 2798, 1[i0[], 2[3[i1[], i2[]], i3[]]]: 2782, 1[i0[], 2[i3[], 3[i2[], i1[]]]]: 2741, 1[3[i0[], i1[]], 2[i2[], i3[]]]: 2754, 1[2[i0[], 3[i3[], i1[]]], i2[]]: 2824, 1[2[i0[], i2[]], 3[i3[], i1[]]]: 2792, 1[i0[], 2[i2[], 3[i3[], i1[]]]]: 2865, 1[2[3[i0[], i1[]], i2[]], i3[]]: 2913, 1[3[i0[], i3[]], 2[i1[], i2[]]]: 2849, 1[2[i0[], 3[i3[], i2[]]], i1[]]: 2661, 1[i0[], 2[i1[], 3[i2[], i3[]]]]: 2792, 1[2[3[i0[], i2[]], i1[]], i3[]]: 2740, 1[i0[], 2[3[i2[], i3[]], i1[]]]: 2780, 1[2[i0[], i3[]], 3[i1[], i2[]]]: 2689, 1[2[3[i0[], i3[]], i2[]], i1[]]: 2764, 1[2[3[i0[], i2[]], i3[]], i1[]]: 2769, 1[i0[], 2[3[i2[], i1[]], i3[]]]: 2731, 1[2[3[i0[], i3[]], i1[]], i2[]]: 2826}\n"}︡ ︠9874a71e-dd7a-44a9-aff4-1def37046de0︠ d=CountsDict([forgetLeafLabels(t) for t in ts]) d ︡6fbbab3a-e0ff-4579-9132-37793d142aa0︡{"stdout":"{1[0[], 2[3[0[], 0[]], 0[]]]: 16634, 1[0[], 2[0[], 3[0[], 0[]]]]: 16719, 1[2[3[0[], 0[]], 0[]], 0[]]: 16782, 1[3[0[], 0[]], 2[0[], 0[]]]: 16769, 1[2[0[], 0[]], 3[0[], 0[]]]: 16658, 1[2[0[], 3[0[], 0[]]], 0[]]: 16438}\n"}︡ ︠1797864f-906e-42c5-a5b3-6ce6e04b9d1e︠ d=CountsDict([forgetAllLabels(t) for t in ts]) d ︡61485a83-f55c-4a05-9340-ddfa6ad2e2ab︡{"stdout":"{[[[], []], [[], []]]: 33427, [[[[], []], []], []]: 16782, [[], [[[], []], []]]: 16634, [[], [[], [[], []]]]: 16719, [[[], [[], []]], []]: 16438}\n"}︡ ︠b36d4b84-7105-4c94-a178-2300bffc7738︠ d=CountsDict([justTree(t) for t in ts]) d ︡14552da7-32d2-44c8-bcb0-fa956a124578︡{"stdout":"{Graph on 7 vertices: 66573, Graph on 7 vertices: 33427}\n"}︡ ︠07c1ab4f-f0c3-44ce-ab03-98d0cedd17be︠ graphs.CompleteGraph(5).to_directed().show() ︡f978c906-ff72-4040-8c1e-c7b1870e03d7︡{"once":false,"file":{"show":true,"uuid":"ba831715-60b8-4756-b650-6b47685994a5","filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute1-us/16258/tmp_oI3Kv6.svg"}}︡{"html":""}︡ ︠fdee4dfb-47b3-49bd-9f0f-431eb036858b︠ ts=[transmissionProcessTC(graphs.CompleteGraph(5).to_directed(),0) for _ in range(100000)] d=CountsDict([forgetLeafLabels(t) for t in ts]) print len(d) d ︡d281780e-0491-45a1-8f6c-e31182b4a992︡{"stdout":"24\n"}︡{"stdout":"{1[3[0[], 0[]], 2[4[0[], 0[]], 0[]]]: 4283, 1[0[], 2[3[0[], 0[]], 4[0[], 0[]]]]: 4119, 1[2[4[0[], 0[]], 0[]], 3[0[], 0[]]]: 4232, 1[3[0[], 4[0[], 0[]]], 2[0[], 0[]]]: 4124, 1[0[], 2[0[], 3[4[0[], 0[]], 0[]]]]: 4110, 1[3[4[0[], 0[]], 0[]], 2[0[], 0[]]]: 4164, 1[2[0[], 4[0[], 0[]]], 3[0[], 0[]]]: 4244, 1[4[0[], 0[]], 2[3[0[], 0[]], 0[]]]: 4195, 1[2[0[], 3[0[], 4[0[], 0[]]]], 0[]]: 4158, 1[2[0[], 0[]], 3[4[0[], 0[]], 0[]]]: 3997, 1[2[3[0[], 0[]], 4[0[], 0[]]], 0[]]: 4094, 1[2[4[0[], 0[]], 3[0[], 0[]]], 0[]]: 4200, 1[0[], 2[4[0[], 0[]], 3[0[], 0[]]]]: 4186, 1[2[0[], 3[0[], 0[]]], 4[0[], 0[]]]: 4170, 1[2[3[0[], 4[0[], 0[]]], 0[]], 0[]]: 4205, 1[0[], 2[3[0[], 4[0[], 0[]]], 0[]]]: 4165, 1[0[], 2[3[4[0[], 0[]], 0[]], 0[]]]: 4109, 1[2[3[0[], 0[]], 0[]], 4[0[], 0[]]]: 4188, 1[0[], 2[0[], 3[0[], 4[0[], 0[]]]]]: 4179, 1[2[3[4[0[], 0[]], 0[]], 0[]], 0[]]: 4222, 1[4[0[], 0[]], 2[0[], 3[0[], 0[]]]]: 4081, 1[2[0[], 3[4[0[], 0[]], 0[]]], 0[]]: 4212, 1[2[0[], 0[]], 3[0[], 4[0[], 0[]]]]: 4186, 1[3[0[], 0[]], 2[0[], 4[0[], 0[]]]]: 4177}\n"}︡ ︠11c0574f-0173-4638-9035-d27509ae6d3a︠ d=CountsDict([forgetAllLabels(t) for t in ts]) print len(d) d ︡d9e85b92-083d-4da5-94a3-0f8fccafe65c︡{"stdout":"14\n"}︡{"stdout":"{[[[[], [[], []]], []], []]: 4205, [[], [[[], []], [[], []]]]: 8305, [[[], []], [[], [[], []]]]: 12444, [[[], [[[], []], []]], []]: 4212, [[[[], []], [[], []]], []]: 8294, [[[], [[], []]], [[], []]]: 12538, [[], [[[[], []], []], []]]: 4109, [[[[], []], []], [[], []]]: 12584, [[], [[], [[], [[], []]]]]: 4179, [[[], []], [[[], []], []]]: 12475, [[[[[], []], []], []], []]: 4222, [[[], [[], [[], []]]], []]: 4158, [[], [[[], [[], []]], []]]: 4165, [[], [[], [[[], []], []]]]: 4110}\n"}︡ ︠c718e0b0-b290-431a-8542-d6a30e38bf91︠ d=CountsDict([justTree(t) for t in ts]) len(d) d ︡87dbfcd8-2b5c-4121-9070-1d4fb233f990︡{"stdout":"3\n"}︡{"stdout":"{Graph on 9 vertices: 33360, Graph on 9 vertices: 50041, Graph on 9 vertices: 16599}\n"}︡ ︠52a09746-cbdf-42ed-aba2-a5b9a712f61bi︠ %htmlLet the SICN be the path graph. We get the distribution of transmission trees at various resolutions next.
︡d0470eb6-021c-4946-958d-cc01758d94ee︡{"html":"\n Let the SICN be the path graph. We get the distribution of transmission trees at various resolutions next.\n
\n"}︡ ︠fe6f2726-2a3d-41fe-8fd1-3a2c35d64797︠ digraphs.Path(4).show() ︡a35e29a7-d81a-4aa8-98ae-fb9bf9a019ea︡{"once":false,"file":{"show":true,"uuid":"b33932f3-44e3-4ac2-a194-276d1f77baea","filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute1-us/16258/tmp_b_jHtq.svg"}}︡{"html":""}︡ ︠54b7841a-1e77-46c7-9f8c-6b2889e06d84︠ # demo transmissionProcessTC(digraphs.Path(4),0) ︡695991a9-7154-47e0-ab47-f06fa07a6d1a︡{"stdout":"step z = 1\n 1_\n / \\\ni0 i1\n--------------------\nstep z = 2\n _1__\n / \\\ni0 2_\n / \\\n i1 i2\n--------------------\nstep z = 3\n __1__\n / \\\ni0 _2__\n / \\\n i1 3_\n / \\\n i2 i3\n--------------------\n1[i0[], 2[i1[], 3[i2[], i3[]]]]\n"}︡ ︠77b29186-6718-4a9e-8007-8db5f5d077f4︠ ts=[transmissionProcessTC(digraphs.Path(4),0) for _ in range(1000)] d=CountsDict(ts) print len(d) d ︡67e31590-3fb3-42b3-9253-afeab359ac7c︡{"stdout":"1\n"}︡{"stdout":"{1[i0[], 2[i1[], 3[i2[], i3[]]]]: 1000}\n"}︡ ︠0af39fec-0464-4fd8-9ac8-a5a1c493bfe4︠ d=CountsDict([forgetAllLabels(t) for t in ts]) print len(d) d ︡f21e322e-950c-41df-b31c-58ad93182a46︡{"stdout":"1\n"}︡{"stdout":"{[[], [[], [[], []]]]: 1000}\n"}︡ ︠05a460a1-84e1-4a5b-891b-34921ce9f8ed︠ d=CountsDict([justTree(t) for t in ts]) len(d) d ︡8de78073-8df3-499c-8828-4280522e5ff8︡{"stdout":"1\n"}︡{"stdout":"{Graph on 7 vertices: 1000}\n"}︡ ︠5ac817a1-6ec0-47e8-a7af-2d4f0f56c53di︠ %htmlLet the SICN be the star network. We get the distribution of transmission trees at various resolutions next.
︡88f07a1c-7b4f-4027-906d-a5bc6a01035d︡{"html":"\n Let the SICN be the star network. We get the distribution of transmission trees at various resolutions next.\n
\n"}︡ ︠7755fe00-bf23-4069-810a-562dc97a9b82︠ g=graphs.StarGraph(3).to_directed() g.show() ︡220129e3-3831-4bab-923d-4ec37b7f84ea︡{"once":false,"file":{"show":true,"uuid":"9c583efc-8d35-47f8-a9b4-fe4d2970c0e9","filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute1-us/16258/tmp_cX0r_A.svg"}}︡{"html":""}︡ ︠49de74ad-2572-418f-a750-2b9decb4a6aa︠ # demo transmissionProcessTC(graphs.StarGraph(3).to_directed(),0) ︡25aa7f13-f1f9-4272-a4c8-bc2364514387︡{"stdout":"step z = 1\n 1_\n / \\\ni0 i3\n--------------------\nstep z = 2\n __1__\n / \\\n 2_ i3\n / \\ \ni0 i2 \n--------------------\nstep z = 3\n __1___\n / \\\n __2__ i3\n / \\ \n 3_ i2 \n / \\ \ni0 i1 \n--------------------\n1[2[3[i0[], i1[]], i2[]], i3[]]\n"}︡ ︠486c9871-6dca-4011-96dd-fdf38dbf5900︠ ts=[transmissionProcessTC(graphs.StarGraph(3).to_directed(),0) for _ in range(10000)] d=CountsDict(ts) print len(d) d ︡633bab1e-4f10-425e-834c-d891b5767c62︡{"stdout":"6\n"}︡{"stdout":"{1[2[3[i0[], i2[]], i1[]], i3[]]: 1633, 1[2[3[i0[], i1[]], i2[]], i3[]]: 1645, 1[2[3[i0[], i3[]], i2[]], i1[]]: 1692, 1[2[3[i0[], i2[]], i3[]], i1[]]: 1682, 1[2[3[i0[], i1[]], i3[]], i2[]]: 1674, 1[2[3[i0[], i3[]], i1[]], i2[]]: 1674}\n"}︡ ︠4d1d334a-7cf5-47c5-8b97-8988a8e3d1b5︠ d=CountsDict([forgetLeafLabels(t) for t in ts]) print len(d) d ︡6921e244-d762-48b9-83fd-d42fd63f7815︡{"stdout":"1\n"}︡{"stdout":"{1[2[3[0[], 0[]], 0[]], 0[]]: 10000}\n"}︡ ︠6557ab59-c63b-417d-be7a-66b68e22b434︠ d=CountsDict([forgetAllLabels(t) for t in ts]) print len(d) d ︡70d11285-911c-4e88-9eb9-dc915a449acc︡{"stdout":"1\n"}︡{"stdout":"{[[[[], []], []], []]: 10000}\n"}︡ ︠96d991d8-ff79-4509-88a0-fea6aa33169e︠ d=CountsDict([justTree(t) for t in ts]) len(d) d ︡4a922824-33f7-416c-b428-984d8bf4c3b8︡{"stdout":"1\n"}︡{"stdout":"{Graph on 7 vertices: 10000}\n"}︡ ︠e2ee1627-bfc5-4504-b3cd-7421eeb067b0i︠ %md # Likelihood Functions & Sufficient Statistics ## based on tree topology without labels or branch-lengths Note one can include branch-lengths into the likelihood easily via the exponential rates for the continuous time Markov processes that the discrete jump chains studied here. This will be necessary to conduct statistical inference for applied epidemiological models. ︡76f73d2b-df57-4ee3-8d05-3cd8606a443b︡{"done":true,"md":"# Likelihood Functions & Sufficient Statistics\n## based on tree topology without labels or branch-lengths\nNote one can include branch-lengths into the likelihood easily via the exponential rates for the continuous time Markov processes that the discrete jump chains studied here.\nThis will be necessary to conduct statistical inference for applied epidemiological models."} ︠58213fa6-6edf-4329-adfa-2d29bae92fdfs︠ def splitsSequence(T): '''return a list of tuples (left,right) split sizes at each split node''' l = [] LabelledBinaryTree(T).post_order_traversal(lambda node: l.append((node[0].node_number(),node[1].node_number()))) return l def prob_RPT(T,a,b): '''probability of ranked planar tree T under beta-splitting model a,b>-1, where (a+1,b+1) are the parameters of the beta distribution''' return prod(map(lambda x: beta(x[0]+a+1,x[1]+b+1)/beta(a+1,b+1), splitsSequence(T))) def negLogLkl_SplitPairCounts(spc,a,b): '''-log likelihood of multiple independent ranked planar trees through their sufficient statistics of the frequence of split-pair counts spc= [(nL_i,nR_i,c_i): i=1,..,K] under beta-splitting model a,b>-1, where (a+1,b+1) are the parameters of the beta distribution -- This implements first Equation in Thm 1 involving beta functions''' return -RR(sum(map(lambda x: x[2]*log(1.0*beta(x[0]+a+1,x[1]+b+1)/beta(a+1,b+1)), spc))) def splitPairsCounts(TS): '''list of the frequency of all distinct split-pairs, i.e. (# of left splits, # right splits) beloe each internal node in each transmission tree in the list TS of transmission trees''' splitPairCounts=sorted(CountsDict(flatten([splitsSequence(t) for t in TS],max_level=1)).items()) return [(x[0][0],x[0][1],x[1]) for x in splitPairCounts] def splitPairsCountsDict(TS): '''dictionary of the frequency of all distinct split-pairs, i.e. (# of left splits, # right splits) beloe each internal node in each transmission tree in the list TS of transmission trees''' #splitPairCounts=sorted(CountsDict(flatten([splitsSequence(t) for t in TS],max_level=1))) sD = CountsDict(flatten([splitsSequence(t) for t in TS],max_level=1)) return sD def logLklOfASplitPair(a,b,nL,nR): '''beta(nL+a+1,nR+b+1)/beta(a+1,b+1) without beta functions via Eqn 2 in Thm 1''' A1=sum([log((b+j)/(b+j+a)) for j in range(nR+1)]) A2=sum([log((a+i)/(a+i+b+nR+1)) for i in range(nL+1)]) A3=log(b*a/((a+b)*(a+b+1))) return A1+A2-A3 def negLogLkl_SplitPairCounts2(spc,a,b): '''-log likelihood of multiple independent ranked planar trees through their sufficient statistics of the frequency of split-pair counts spc= [(nL_i,nR_i,c_i): i=1,..,K] under beta-splitting model a,b>-1, where (a+1,b+1) are the parameters of the beta distribution -- This implements second Equation in Thm 1 without beta functions''' return -(sum(map(lambda x: x[2]*logLklOfASplitPair(a,b,x[0],x[1]), spc))) def LklOfASplitPair(a,b,nL,nR): '''beta(nL+a+1,nR+b+1)/beta(a+1,b+1) without beta functions via Eqn 2 in Thm 1''' A1=prod([((b+j)/(b+j+a)) for j in range(nR+1)]) A2=prod([((a+i)/(a+i+b+nR+1)) for i in range(nL+1)]) A3=(b*a/((a+b)*(a+b+1))) return (A1*A2)/A3 def negLogLkl_SplitPairCounts2Prod(spc,a,b): '''- log likelihood of multiple independent ranked planar trees through their sufficient statistics of the frequency of split-pair counts spc= [(nL_i,nR_i,c_i): i=1,..,K] under beta-splitting model a,b>-1, where (a+1,b+1) are the parameters of the beta distribution -- This implements second Equation in Thm 1 without beta functions''' return -(sum(map(lambda x: x[2]*log(LklOfASplitPair(a,b,x[0],x[1])), spc))) ︡e8db63d5-862b-404e-a16e-245c5cce2a0a︡{"done":true}︡ ︠42bddcab-b6e9-48ee-bd8b-a9e4300805aci︠ %md ## Demo of the mle for a complete graph with 50 nodes and 10 sampled trees ︡55a1f64e-a8a9-4795-84d0-d7daa6728777︡{"done":true,"md":"## Demo of the mle for a complete graph with 50 nodes and 10 sampled trees"} ︠e5dad933-edbb-49d3-a234-50959dfd28bd︠ %time # demo of the mle for a complete graph with 50 nodes and 10 sampled trees c_1 = lambda p: p[0]+0.9999999 # constraint for alpha > -1 c_2 = lambda p: p[1]+0.9999999 # constraint for beta > -1 # simulation settings n=50 reps=10 # trial 1: make a list of 10 independent transmission trees on complete graph with 50 nodes ts1=[transmissionProcessTC(graphs.CompleteGraph(n).to_directed(),0) for _ in range(reps)] spc1=splitPairsCounts(ts1) # trial 2: make another list of 10 independent transmission trees on complete graph with 50 nodes ts2=[transmissionProcessTC(graphs.CompleteGraph(n).to_directed(),0) for _ in range(reps)] spc2=splitPairsCounts(ts2) ︡7880a133-e70d-4bd6-a586-e480dc461c31︡{"stdout":"CPU time: 1.94 s, Wall time: 3.89 s\n"}︡{"done":true}︡ ︠b3f97b9b-a7b8-4cad-b02c-ce10df6e6bdf︠ spc1=splitPairsCounts(ts1) ︡0a789876-8fcb-4cfe-8db1-d8aa4861c7b9︡{"done":true}︡ ︠2d114f33-6bdb-4ea1-83fa-466fa6d2f763︠ # slower more explicit method of finding MLE # (don't use this for larger simuations!) def negLkl1(AB): # negative log likelihood function for trees from trial 1 return sum([-log(1.0*prob_RPT(ts1[j],AB[0],AB[1])) for j in range(reps)]) def negLkl2(AB): # negative log likelihood function for trees from trial 1 return sum([-log(1.0*prob_RPT(ts2[j],AB[0],AB[1])) for j in range(reps)]) %time mle=minimize_constrained(negLkl1,[c_1,c_2],[0.0,0.0],disp=0) #MLE for trial 1 print [n,reps,mle] %time mle=minimize_constrained(negLkl2,[c_1,c_2],[0.0,0.0],disp=0) #MLE for trial 2 print [n,reps,mle] ︡a5f99cfa-81b0-4d15-b223-baaf62f0d637︡{"stdout":"CPU time: 49.64 s, Wall time: 100.20 s"}︡{"stdout":"\n"}︡{"stdout":"[50, 10, (-0.06665985295552748, -0.05038096402831602)]\n"}︡{"stdout":"CPU time: 41.69 s, Wall time: 84.02 s"}︡{"stdout":"\n"}︡{"stdout":"[50, 10, (0.004389330784528777, -0.04321792852343554)]\n"}︡{"done":true}︡ ︠5deb1d1d-0de1-4d70-9a85-83e0097dacfb︠ # faster MLE computation using sufficient statistics of split-pair frequencies # (don't use this for larger simuations!) def negLkl1(AB): return negLogLkl_SplitPairCounts(spc1,AB[0],AB[1]) def negLkl2(AB): return negLogLkl_SplitPairCounts(spc2,AB[0],AB[1]) %time mle=minimize_constrained(negLkl1,[c_1,c_2],[0.0,0.0],disp=0) #MLE for trial 1 print [n,reps,mle] %time mle=minimize_constrained(negLkl2,[c_1,c_2],[0.0,0.0],disp=0) #MLE for trial 2 print [n,reps,mle] ︡e24f8cf7-c023-4556-a5dd-d76b82d34621︡{"stdout":"CPU time: 9.43 s, Wall time: 19.07 s"}︡{"stdout":"\n"}︡{"stdout":"[50, 10, (-0.06665985290623813, -0.05038096406044273)]\n"}︡{"stdout":"CPU time: 7.18 s, Wall time: 14.59 s"}︡{"stdout":"\n"}︡{"stdout":"[50, 10, (0.004740228015888741, -0.04293877679246926)]\n"}︡{"done":true}︡ ︠52ebcd0a-6575-4381-a0fe-8638fcc040c7︠ # fastest and numerically most robust MLE computation using sufficient statistics # USE this for larger simulations def negLkl1(AB): return negLogLkl_SplitPairCounts2(spc1,AB[0],AB[1]) def negLkl2(AB): return negLogLkl_SplitPairCounts2(spc2,AB[0],AB[1]) %time mle=minimize_constrained(negLkl1,[c_1,c_2],[0.0,0.0],disp=0) #MLE for trial 1 print [n,reps,mle] %time mle=minimize_constrained(negLkl2,[c_1,c_2],[0.0,0.0],disp=0) #MLE for trial 2 print [n,reps,mle] ︡42cb8f35-385e-4886-9605-44854ecd0981︡{"stdout":"CPU time: 4.73 s, Wall time: 9.52 s"}︡{"stdout":"\n"}︡{"stdout":"[50, 10, (-0.06644703305884028, -0.05022885910313697)]\n"}︡{"stdout":"CPU time: 9.34 s, Wall time: 18.87 s"}︡{"stdout":"\n"}︡{"stdout":"[50, 10, (0.004665943054046424, -0.043007797370934346)]\n"}︡{"done":true}︡ ︠49653165-aba2-416f-8423-04107d31d200i︠ %md ## Sufficient Statistics of the parameters specifying the distribution for Beta-splitting Transmission Tree topologies ︡4e7722b9-75d4-42e6-8f43-c8a82be05d48︡{"done":true,"md":"## Sufficient Statistics of the parameters specifying the distribution for Beta-splitting Transmission Tree topologies"} ︠f907ddca-1224-40d8-b4be-206db7f66783︠ len(ts1),len(ts2) # number of trees in each trial ︡4866eb71-3a3b-42ee-89b2-a69b3d36ad54︡{"stdout":"(10, 10)\n"}︡{"done":true}︡ ︠3817f021-451e-4b43-82ac-9db453ea0952︠ def denMatOfsplitPairsCounts(ts,n,EMF): '''dense matrix of split pair counts of transmission trees in the list ts normalise to sum to 1 if EMF = Empirical Mass Function = True ''' spcD=splitPairsCountsDict(ts); spcD[(0,n)]=0; spcD[(n,0)]=0; ss=matrix(spcD).dense_matrix() if EMF: ss=ss/sum(sum(ss)) return ss ︡e7291345-044e-4fbb-8a2e-80bca33944aa︡{"done":true}︡ ︠5d9d4589-dc71-4324-9054-da0e16c4ec42i︠ %md Let's save these plots ︡27b84dde-616d-4673-9e43-7425110a432c︡{"done":true,"md":"Let's save these plots"} ︠bb364dc1-f884-41db-8ec7-cac6ac98c807︠ mp1=matrix_plot(np.log1p(denMatOfsplitPairsCounts(ts1,n,0)) , cmap='summer', colorbar=True) mp2=matrix_plot(np.log1p(denMatOfsplitPairsCounts(ts2,n,0)) , cmap='summer', colorbar=True) ︡3f17e3c9-1690-4599-a042-b6856b85fd50︡{"done":true}︡ ︠a50bdc58-ffab-41de-baa0-99675e0ed40c︠ mp2.show() ︡6c32f69c-d376-485c-9c76-f2eabe1ec6a6︡{"file":{"filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute5-us/21948/tmp_4i8PGp.svg","show":true,"text":null,"uuid":"a95340d8-a412-466b-97f2-470bd32c8959"},"once":false}︡{"html":""}︡{"done":true}︡ ︠5dfbf597-ad91-4f24-9963-43fb96c26eaa︠ %sh pwd ︡6b38e19d-48b1-40b3-9c3d-81001de147fb︡{"stdout":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c\n"}︡{"done":true}︡ ︠fb1f0962-e916-4b96-b48e-475a7dfde9d7︠ %sh #mkdir /projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/figures && ls ls ︡1bf3af2b-b336-416b-86da-b79ef21ea742︡{"stdout":"2015-10-25-165503.sagews\n2016-05-30-transmissionProcess.sagews\nfigures\n"}︡{"done":true}︡ ︠13c2e970-88f7-4843-bbf6-af906aee3837︠ mp1.save('./figures/SuffStatsCompleteGraph_n50_reps10_mp1.png') mp1.save('./figures/SuffStatsCompleteGraph_n50_reps10_mp1.eps') mp1.save('./figures/SuffStatsCompleteGraph_n50_reps10_mp1.pdf') mp2.save('./figures/SuffStatsCompleteGraph_n50_reps10_mp2.png') mp2.save('./figures/SuffStatsCompleteGraph_n50_reps10_mp2.eps') mp2.save('./figures/SuffStatsCompleteGraph_n50_reps10_mp2.pdf') ︡84b7ea86-263d-4fdc-bffd-afc9152672a0︡{"done":true}︡ ︠043267d5-e5f7-4743-8339-4c5a31362ccdi︠ %md Let's do a more exhaustive simulation next ︡447aae22-a625-4b1e-aa9b-240b5c4053ba︡{"done":true,"md":"Let's do a more exhaustive simulation next"} ︠ca166c69-4fd6-4f56-8187-ff5e50715330︠ %time # demo of the mle for a complete graph with 50 nodes and 10 sampled trees c_1 = lambda p: p[0]+0.9999999 # constraint for alpha > -1 c_2 = lambda p: p[1]+0.9999999 # constraint for beta > -1 # simulation settings n=50 reps=1000 # trial 1: make a list of 10 independent transmission trees on complete graph with 50 nodes ts3=[transmissionProcessTC(graphs.CompleteGraph(n).to_directed(),0) for _ in range(reps)] spc3=splitPairsCounts(ts3) ︡a8292275-9bc3-4e38-818b-10cc21c6fb88︡{"stdout":"CPU time: 74.06 s, Wall time: 150.04 s"}︡{"stdout":"\n"}︡{"done":true}︡ ︠a26493d7-9e2f-4b57-816f-173402ba087c︠ spc3=splitPairsCounts(ts3) ︡7939aff6-6df9-4b3d-ae69-f9ea08c6806b︡{"done":true}︡ ︠6d7e6229-f1ec-4d6e-aa2e-bd591b2b9ff8︠ def negLkl3(AB): return negLogLkl_SplitPairCounts2(spc3,AB[0],AB[1]) %time mle=minimize_constrained(negLkl3,[c_1,c_2],[0.0,0.0],disp=0) #MLE for trial 1 print [n,reps,mle] ︡fd4c9a69-efbc-4643-918d-33445c5260ca︡{"stdout":"CPU time: 76.10 s, Wall time: 153.64 s"}︡{"stdout":"\n"}︡{"stdout":"[50, 1000, (0.027896177617489963, 0.032449437588353475)]\n"}︡{"done":true}︡ ︠6b79916e-80af-473f-9fe0-64def262da41i︠ %md Let's plot the sufficient stats for this longer run. ︡27137afc-0861-4a41-9794-9ffe5e3f11a4︡{"hide":"input"}︡{"md":"Let's plot the sufficient stats for this longer run.\n"}︡{"done":true}︡ ︠a5718b16-3f31-4885-baf5-212969c0c6cb︠ #mp3=matrix_plot(np.log1p(denMatOfsplitPairsCounts(ts3,n,1)*np.double(1000)), norm='value', cmap='summer') # a log(1+prob*1000) plot of the mp3=matrix_plot(np.log1p(denMatOfsplitPairsCounts(ts3,n,0)*np.double(1)) , cmap='summer', colorbar=True)#, norm='value', cmap='summer', colorbar=True) # a log(1+freq) plot of the ︡b3378e3b-0b8b-44ec-bca7-2f32d8c5bc24︡{"done":true}︡ ︠1dabc8cc-2bf3-44f3-ba4f-6940ff9da8e0︠ mp3.show() ︡0958bcbb-b6d4-483d-aded-db672b755c02︡{"file":{"filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute5-us/21948/tmp_oBy9PT.svg","show":true,"text":null,"uuid":"e24ba1ec-6bfc-42ca-8fbf-b5cce31eb81e"},"once":false}︡{"html":""}︡{"done":true}︡ ︠912f5465-3be1-4f23-a4b9-02c2aa8c4476︠ mp3.save('./figures/SuffStatsCompleteGraph_n50_reps1000_mp3.png') mp3.save('./figures/SuffStatsCompleteGraph_n50_reps1000_mp3.eps') mp3.save('./figures/SuffStatsCompleteGraph_n50_reps1000_mp3.pdf') ︡a51abd3f-4c8c-4923-99f6-0b43d3fce09c︡{"done":true}︡ ︠6bf7dfeb-0e67-4be5-bde9-42760d510b84︠ %sh ls figures/ ︡651fd3e0-5637-40e4-9e2d-494bbff06fc7︡{"stdout":"SuffStatsCompleteGraph_n50_reps10_mp1.eps\nSuffStatsCompleteGraph_n50_reps10_mp1.pdf\nSuffStatsCompleteGraph_n50_reps10_mp1.png\nSuffStatsCompleteGraph_n50_reps10_mp2.eps\nSuffStatsCompleteGraph_n50_reps10_mp2.pdf\nSuffStatsCompleteGraph_n50_reps10_mp2.png\n"}︡{"done":true}︡ ︠82f2922e-95eb-4844-a96d-22025bfa5d12i︠ %md # Some Deterministic Contact Networks ︡57b882d6-9449-44fd-aa8f-81c4a5a0d6ec︡{"done":true,"md":"# Some Deterministic Contact Networks"} ︠2f231a43-4096-4c7d-9f44-34f9bb1d7502i︠ %md ## Bidirectional Circular Network ︡1e63ccdc-9b38-4a47-99fc-b9e3d5689810︡{"done":true,"md":"## Bidirectional Circular Network"} ︠152bdc65-306d-4278-b2ec-899af116ac7es︠ g=digraphs.Circulant(5,[1,-1]) g.show() ︡01a46ff3-5eac-4c73-9bd4-91ad12e1f453︡{"file":{"filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute5-us/23053/tmp_AfMclJ.svg","show":true,"text":null,"uuid":"f0adbdec-ff64-49c5-936a-a80b040874d5"},"once":false}︡{"html":""}︡{"done":true}︡ ︠314cc338-f44b-4cfb-ae56-8e128f422cfci︠ %md To Make Latex tikz-graph see [http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html](http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html). ︡734341e0-7b0d-4dd7-b97c-5d1d9feebbfc︡{"done":true,"md":"To Make Latex tikz-graph see [http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html](http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html)."} ︠241d1d8a-911b-4ec1-81c5-c1836cf57810s︠ H = digraphs.Circulant(6,[1,-1]) H.set_latex_options( graphic_size=(5,5), vertex_size=0.2, edge_thickness=0.04, edge_color='green', vertex_color='green', vertex_label_color='red', format='tkz_graph', #tkz_style='Art', #layout='acyclic' ) latex(H) ︡bd2fcd58-e7b4-4ba9-ba3e-279ce15f3631︡{"stdout":"\\begin{tikzpicture}\n%\n\\useasboundingbox (0,0) rectangle (5.0cm,5.0cm);\n%\n\\definecolor{cv0}{rgb}{0.0,0.502,0.0}\n\\definecolor{cfv0}{rgb}{1.0,1.0,1.0}\n\\definecolor{clv0}{rgb}{1.0,0.0,0.0}\n\\definecolor{cv1}{rgb}{0.0,0.502,0.0}\n\\definecolor{cfv1}{rgb}{1.0,1.0,1.0}\n\\definecolor{clv1}{rgb}{1.0,0.0,0.0}\n\\definecolor{cv2}{rgb}{0.0,0.502,0.0}\n\\definecolor{cfv2}{rgb}{1.0,1.0,1.0}\n\\definecolor{clv2}{rgb}{1.0,0.0,0.0}\n\\definecolor{cv3}{rgb}{0.0,0.502,0.0}\n\\definecolor{cfv3}{rgb}{1.0,1.0,1.0}\n\\definecolor{clv3}{rgb}{1.0,0.0,0.0}\n\\definecolor{cv4}{rgb}{0.0,0.502,0.0}\n\\definecolor{cfv4}{rgb}{1.0,1.0,1.0}\n\\definecolor{clv4}{rgb}{1.0,0.0,0.0}\n\\definecolor{cv5}{rgb}{0.0,0.502,0.0}\n\\definecolor{cfv5}{rgb}{1.0,1.0,1.0}\n\\definecolor{clv5}{rgb}{1.0,0.0,0.0}\n\\definecolor{cv0v1}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv0v5}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv1v0}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv1v2}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv2v1}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv2v3}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv3v2}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv3v4}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv4v3}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv4v5}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv5v0}{rgb}{0.0,0.502,0.0}\n\\definecolor{cv5v4}{rgb}{0.0,0.502,0.0}\n%\n\\Vertex[style={minimum size=0.2cm,draw=cv0,fill=cfv0,text=clv0,shape=circle},LabelOut=false,L=\\hbox{$0$},x=5.0cm,y=2.5cm]{v0}\n\\Vertex[style={minimum size=0.2cm,draw=cv1,fill=cfv1,text=clv1,shape=circle},LabelOut=false,L=\\hbox{$1$},x=3.75cm,y=5.0cm]{v1}\n\\Vertex[style={minimum size=0.2cm,draw=cv2,fill=cfv2,text=clv2,shape=circle},LabelOut=false,L=\\hbox{$2$},x=1.25cm,y=5.0cm]{v2}\n\\Vertex[style={minimum size=0.2cm,draw=cv3,fill=cfv3,text=clv3,shape=circle},LabelOut=false,L=\\hbox{$3$},x=0.0cm,y=2.5cm]{v3}\n\\Vertex[style={minimum size=0.2cm,draw=cv4,fill=cfv4,text=clv4,shape=circle},LabelOut=false,L=\\hbox{$4$},x=1.25cm,y=0.0cm]{v4}\n\\Vertex[style={minimum size=0.2cm,draw=cv5,fill=cfv5,text=clv5,shape=circle},LabelOut=false,L=\\hbox{$5$},x=3.75cm,y=0.0cm]{v5}\n%\n\\Edge[lw=0.04cm,style={post, bend right,color=cv0v1,},](v0)(v1)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv0v5,},](v0)(v5)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv1v0,},](v1)(v0)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv1v2,},](v1)(v2)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv2v1,},](v2)(v1)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv2v3,},](v2)(v3)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv3v2,},](v3)(v2)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv3v4,},](v3)(v4)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv4v3,},](v4)(v3)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv4v5,},](v4)(v5)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv5v0,},](v5)(v0)\n\\Edge[lw=0.04cm,style={post, bend right,color=cv5v4,},](v5)(v4)\n%\n\\end{tikzpicture}"}︡{"stdout":"\n"}︡{"done":true}︡ ︠32111e48-685b-4f31-b52f-e3f4788b0ff5︠ ts=[transmissionProcessTC(digraphs.Circulant(6,[1,-1]),0) for _ in range(10)] ︡790b9926-198b-42b8-8153-47f3d9f34a43︡{"done":true}︡ ︠1ef3c039-2dc2-4dc5-a113-213c5da59b9e︠ T = ts[3] view(T) ︡30ac5c99-4e33-4813-bc84-13419ba3d850︡{"file":{"filename":"/tmp/tmp23pyy2.png","show":true,"text":null,"uuid":"a3adda8e-b3ce-4f4f-bd17-5a4d2310863b"},"once":false}︡{"html":""}︡{"done":true}︡ ︠69bd1673-c000-4e93-a3b1-c14eff79468d︠ #latex(T) ︡ccaa57e3-63df-4327-b8c3-f6dfb634e201︡{"done":true}︡ ︠5f295fc2-8d91-4153-be18-48ee2032ff90i︠ %md Just checking that the trees are uniformly distributed over 2^{n-1} ORBOR trees := Only-Right-Branching-subtrees-Of-Root-node trees (for peace of mind :). ︡8bff2ba0-0b8a-4a22-bbcf-b929be7fade0︡{"done":true,"md":"Just checking that the trees are uniformly distributed over 2^{n-1} ORBOR trees := Only-Right-Branching-subtrees-Of-Root-node trees (for peace of mind :)."} ︠5e40c3f4-aff5-4116-815b-ea94f5c1f207s︠ ts=[transmissionProcessTC(digraphs.Circulant(3,[1,-1]),0) for _ in range(10000)] d=CountsDict(ts) print len(d) d ︡2a9d242f-6380-4ac6-9e6f-78d7ea6ada88︡{"stdout":"4\n"}︡{"stdout":"{1[2[i0[], i2[]], i1[]]: 2598, 1[2[i0[], i1[]], i2[]]: 2465, 1[i0[], 2[i2[], i1[]]]: 2492, 1[i0[], 2[i1[], i2[]]]: 2445}\n"}︡{"done":true}︡ ︠eb0bcb7e-6330-4a4a-9a1c-869c45604259s︠ ts=[transmissionProcessTC(digraphs.Circulant(4,[1,-1]),0) for _ in range(10000)] d=CountsDict(ts) print len(d) d ︡93c55822-4487-4033-be10-971ae0e724f6︡{"stdout":"8\n"}︡{"stdout":"{1[2[i0[], i1[]], 3[i3[], i2[]]]: 1264, 1[i0[], 2[i1[], 3[i2[], i3[]]]]: 1226, 1[i0[], 2[i3[], 3[i2[], i1[]]]]: 1255, 1[2[i0[], 3[i1[], i2[]]], i3[]]: 1240, 1[3[i0[], i3[]], 2[i1[], i2[]]]: 1264, 1[2[i0[], 3[i3[], i2[]]], i1[]]: 1240, 1[3[i0[], i1[]], 2[i3[], i2[]]]: 1324, 1[2[i0[], i3[]], 3[i1[], i2[]]]: 1187}\n"}︡{"done":true}︡ ︠2dea160d-378e-4c2a-84f5-cd7234d26810︠ for t in d.keys(): ascii_art(t) ︡34f64246-f67a-4521-9149-e52335d526f7︡{"stdout":" __1____\n / / \n 2__ 3__\n / / / / \ni0 i1 i3 i2\n __1___\n / / \ni0 _2___\n / / \n i1 3__\n / / \n i2 i3\n __1___\n / / \ni0 _2___\n / / \n i3 3__\n / / \n i2 i1\n ___1____\n / / \n _2___ i3\n / / \ni0 3__ \n / / \n i1 i2 \n __1____\n / / \n 3__ 2__\n / / / / \ni0 i3 i1 i2\n ___1____\n / / \n _2___ i1\n / / \ni0 3__ \n / / \n i3 i2 \n __1____\n / / \n 3__ 2__\n / / / / \ni0 i1 i3 i2\n __1____\n / / \n 2__ 3__\n / / / / \ni0 i3 i1 i2\n"}︡{"done":true}︡ ︠a9a85d89-9051-4ee5-927a-f5fd62d9048di︠ %md MLE of transmission trees from the bidirectional circular network. ︡415464dd-5e73-495e-88a6-1e60d38b4981︡{"done":true,"md":"MLE of transmission trees from the bidirectional circular network."} ︠f6380974-6e03-45fc-b74d-a842b2326635︠ %time # demo of the mle for a BalancedTree(2,k) graph with k=4,5,6 or n=31,63,127 nodes and 1 or 10 sampled trees c_1 = lambda p: p[0]+0.9999999 # constraint for alpha > -1 c_2 = lambda p: p[1]+0.9999999 # constraint for beta > -1 mles=[] n=50 reps=1 # all reps have the same unlabelled transmission tree topology trials=5 for i in range(trials): ts=[transmissionProcessTC(digraphs.Circulant(n,[1,-1]),0) for _ in range(reps)] spc=splitPairsCounts(ts) def negLkl(AB): return negLogLkl_SplitPairCounts2(spc,AB[0],AB[1]) mle=minimize_constrained(negLkl,[c_1,c_2],[0.0,0.0],disp=0) x = [n,reps, i, mle, negLkl(mle)] print x mles.append(x) #mles ︡0b4c74a8-3aa4-41e4-88b4-05f4cc4104db︡{"stdout":"[50, 1, 0, (-0.9879913505605317, 1.4153032756705888), 40.135686299805656]"}︡{"stdout":"\n[50, 1, 1, (-0.9874402261518601, 1.56295126382436), 40.71467705498834]"}︡{"stdout":"\n[50, 1, 2, (-0.9878012505019622, 1.5625111252046917), 40.04358386390558]"}︡{"stdout":"\n[50, 1, 3, (-0.9878790181095873, 1.4641878740256553), 40.60995308631155]"}︡{"stdout":"\n[50, 1, 4, (-0.9874316013790505, 1.5708318403090265), 40.781443512465984]"}︡{"stdout":"\n"}︡{"stdout":"CPU time: 33.09 s, Wall time: 33.79 s\n"}︡{"done":true}︡ ︠7173eaf6-81f4-46a6-91a2-151829a50552︠ print (n, reps, trials, mean([x[3][0] for x in mles]), std([x[3][0] for x in mles]), mean([x[3][1] for x in mles]), std([x[3][1] for x in mles])) ︡87c24344-92b0-4345-9e5e-d697bd9ff3a9︡{"stdout":"(50, 1, 5, -0.9880406098481203, 0.0006220653042760461, 1.4583751679644803, 0.15345032109817106)\n"}︡{"done":true}︡ ︠ae0ed783-0002-4cce-95d4-d186818f65e7︠ # (50, 1, 5, -0.9880406098481203, 0.0006220653042760461, 1.4583751679644803, 0.15345032109817106) # (50, 100, 5, -0.9878920259080713, 3.460969702959854e-05, 1.5188788157835482, 0.006721406185256709) ︡8f0fe8b4-50bb-451f-99aa-b4201fc8434a︡ ︠ea8cf74a-8157-41d3-a8f6-f669e8c47dcbi︠ %md Let us repeat the same computation over the relative frequency of split-pairs as opposed to the frequency over a larger number of replicate transmission trees. ︡6fc10327-5d89-47be-a5eb-5b8d5d0daadc︡{"done":true,"md":"Let us repeat the same computation over the relative frequency of split-pairs as opposed to the frequency over a larger number of replicate transmission trees."} ︠38426c5e-c04a-4e63-ade5-9e0b4f567ff2︠ %time # demo of the mle for a BalancedTree(2,k) graph with k=4,5,6 or n=31,63,127 nodes and 1 or 10 sampled trees c_1 = lambda p: p[0]+0.9999999 # constraint for alpha > -1 c_2 = lambda p: p[1]+0.9999999 # constraint for beta > -1 mles=[] n=50 reps=100 # all reps have the same unlabelled transmission tree topology trials=5 for ti in range(trials): ts=[transmissionProcessTC(digraphs.Circulant(n,[1,-1]),0) for _ in range(reps)] spc=splitPairsCounts(ts) s=sum([x[2] for x in spc]) spcRf=[] for i in range(len(spc)): spcRf.append((spc[i][0],spc[i][1],spc[i][2]/s)) def negLkl(AB): return negLogLkl_SplitPairCounts2(spcRf,AB[0],AB[1]) mle=minimize_constrained(negLkl,[c_1,c_2],[0.0,0.0],disp=0) x = [n,reps, i, mle, negLkl(mle)] print x mles.append(x) ︡5bc51ce2-026f-4503-b038-6f4985c490b3︡{"stdout":"[50, 100, 47, (-0.9879013644873202, 1.5266937292679985), 0.8208525765889505]"}︡{"stdout":"\n[50, 100, 52, (-0.9878709556921021, 1.5201429877791235), 0.8183424648045935]"}︡{"stdout":"\n[50, 100, 47, (-0.9880598397493002, 1.4514110195907763), 0.8222490702279951]"}︡{"stdout":"\n[50, 100, 50, (-0.9879350966628759, 1.5017073723940886), 0.8186764955583364]"}︡{"stdout":"\n[50, 100, 47, (-0.98763930925989, 1.517007946626109), 0.8225781541777649]"}︡{"stdout":"\n"}︡{"stdout":"CPU time: 60.00 s, Wall time: 60.83 s\n"}︡{"done":true}︡ ︠3a4c88a1-8ba0-4a2b-b099-acabacfd00b9︠ print (n, reps, trials, mean([x[3][0] for x in mles]), std([x[3][0] for x in mles]), mean([x[3][1] for x in mles]), std([x[3][1] for x in mles])) ︡641cb9f9-7cfb-4392-ade7-264e01633e9d︡{"stdout":"(50, 100, 5, -0.9878813131702978, 0.00015316627305464838, 1.5033926111316191, 0.030470561504457778)\n"}︡{"done":true}︡ ︠8eddd4da-d2d4-41ae-8e49-0653f7295c3e︠ a=-0.98; b=1.5; BetaD = RealDistribution('beta', [a+1, b+1],seed=0) show(BetaD.plot(xmin=0,xmax=1),figsize=[7,2]); ︡d70ac618-1103-4be8-a69e-a2f69aa45601︡{"file":{"filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute5-us/10823/tmp_BfS7E2.svg","show":true,"text":null,"uuid":"d6f37845-9f61-4bcd-815f-959399271979"},"once":false}︡{"html":""}︡{"done":true}︡ ︠3290d61e-a82b-4048-87df-ea748b6d0503︠ ascii_art(ts[4]) # labelled transmission tree ︡76b648af-37de-4788-bc02-a601caf49bee︡{"stdout":" ____1______\n / / \n __2___ i1\n / / \ni0 _3___ \n / / \n i4 4__ \n / / \n i3 i2 \n"}︡{"done":true}︡ ︠5614583e-de1a-4fef-ba7a-a394bfba9687i︠ %md ## Balanced Tree Network ︡8da862be-38aa-4714-9aa6-0a35dda13f48︡{"done":true,"md":"## Balanced Tree Network"} ︠9f9e2306-dec2-4bac-b31e-4fc54558a7fa︠ g=graphs.BalancedTree(3,2).to_directed() g.show() ︡4bcc0af6-ff9d-44c6-9b56-94847f6a3bfe︡{"file":{"filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute5-us/11629/tmp_3jAY62.svg","show":true,"text":null,"uuid":"f0ba7b21-ec34-44bc-a2f1-7838ad05805a"},"once":false}︡{"html":""}︡{"done":true}︡ ︠653b1ed7-b26c-46af-a834-6b5f1611ece4i︠ %md To Make Latex tikz-graph see [http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html](http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html). ︡f56ca2d8-53bb-4d77-af89-290939eeb374︡{"done":true,"md":"To Make Latex tikz-graph see [http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html](http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/graph_latex.html)."} ︠133a3da9-9314-411b-88b4-55550226bef7︠ H = graphs.BalancedTree(2,3).to_directed() #H = graphs.BalancedTree(3,2).to_directed() H.set_latex_options( graphic_size=(5,5), vertex_size=0.2, edge_thickness=0.04, edge_color='green', vertex_color='green', vertex_label_color='red', format='tkz_graph', #tkz_style='Art', layout='acyclic' ) ︡9f670d46-67b3-4a46-a6a4-28cafe844b59︡{"done":true}︡ ︠6215ef10-2bdb-4b23-9140-d0930b125092︠ #view(H) #latex(H) ︡035e3fe1-fe7f-4748-a60b-09bcc5e0ff61︡{"done":true}︡ ︠9d237b20-0a16-4c9b-9618-597566802b96︠ #latex(H) ︡9ec8f40a-87b5-4bab-9836-a6402eea9539︡{"done":true}︡ ︠e56e5b72-10f0-4fbc-8861-99b3420d11eais︠ %md Just checking prob for left-branching d-shark transmission trees (ignoring all labels) under balanced tree network - for peace of mind :) ︡e3a2826f-f78f-40e9-be4a-cae4ed8a5883︡{"hide":"input"}︡{"md":"Just checking prob for left-branching d-shark transmission trees (ignoring all labels) under balanced tree network - for peace of mind :)"}︡{"done":true}︡ ︠4fa0a981-422d-457e-82b2-f8a507b1dedds︠ d,h=3,2 g=graphs.BalancedTree(d,h).to_directed() g.show() ︡e7627e07-0d8e-4dad-b7ca-e3d2f4e97447︡{"file":{"filename":"/projects/58dfa924-55ae-4b6c-9fd4-1cd0ef49eb7c/.sage/temp/compute5-us/23053/tmp_LPZTeS.svg","show":true,"text":null,"uuid":"476b8f41-edbf-4726-8066-d40e0d1c4523"},"once":false}︡{"html":""}︡{"done":true}︡ ︠ad8f5238-c559-4031-a5ee-28bb5016158as︠ d,h=3,3 ts=[transmissionProcessTC(graphs.BalancedTree(d,h).to_directed(),0) for _ in range(10000)] #myDict=CountsDict(ts) #myDict=CountsDict([forgetLeafLabels(t) for t in ts]) myDict=CountsDict([forgetAllLabels(t) for t in ts]) print len(myDict) #myDict ︡82fe2bae-0a71-45b2-81c5-c1365b42f357︡{"stdout":"1\n"}︡{"done":true}︡ ︠c71de5e2-7973-4c7b-8a32-bc8fda3eb615is︠ %md The unlabelled transmission tree is unique. ︡f79a3da2-59e7-4baa-97a4-65666e146ec1︡{"hide":"input"}︡{"md":"The unlabelled transmission tree is unique."}︡{"done":true}︡ ︠3f7b6b7b-05d8-407d-8894-7e73acc89bd1︠ for t in myDict.keys(): # d,h=3,3 "left-branching 3-shark" tree print ascii_art(t) ︡b22e6699-6084-4666-8feb-328eb5329f99︡{"stdout":" _______________________________o________________________________\n / / \n ________________________o__________________________ ________o__________\n / / / / \n ____________o_____________ ________o__________ _______o________ _o__\n / / / / / / / /\no ________o__________ _______o________ _o__ ___o_____ _o__ o__ o\n / / / / / / / / / / / /\n _______o________ _o__ ___o_____ _o__ o__ o o _o__ o__ o o_ o \n / / / / / / / / / / / / / / / / \n ___o_____ _o__ o__ o o _o__ o__ o o_ o o__ o o_ o o o \n / / / / / / / / / / / / / / / / \n o _o__ o__ o o_ o o__ o o_ o o o o_ o o o \n / / / / / / / / / / / / \n o__ o o_ o o o o_ o o o o o \n / / / / / / \n o_ o o o o o \n / / \n o o \n"}︡{"done":true}︡ ︠0b0cf923-a6ef-4f7b-a7e1-ee91196e9aa6︠ # so long so far... could pursue the rank placements on balanced trees ︡83b505fb-c63d-493b-a87c-51f81a00d3bd︡{"done":true}︡ ︠0c43fb10-bdcc-45c8-87e0-78d6efc28788ss︠ graphs.BalancedTree(2,6).to_directed().order() ︡33c6aab8-e0e9-4b89-99fd-bf14cf2fa52d︡{"stdout":"127\n"}︡{"done":true}︡ ︠99f5306d-466c-4cf2-908c-a719976079c7s︠ %time # demo of the mle for a BalancedTree(2,k) graph with k=4,5,6 or n=31,63,127 nodes and 1 or 10 sampled trees c_1 = lambda p: p[0]+0.9999999 # constraint for alpha > -1 c_2 = lambda p: p[1]+0.9999999 # constraint for beta > -1 mles=[] d,h=2,6 n=graphs.BalancedTree(d,h).to_directed().order() reps=1 # all reps have the same unlabelled transmission tree topology trials=2 for i in range(trials): ts=[transmissionProcessTC(graphs.BalancedTree(d,h).to_directed(),0) for _ in range(reps)] spc=splitPairsCounts(ts) def negLkl(AB): return negLogLkl_SplitPairCounts2(spc,AB[0],AB[1]) mle=minimize_constrained(negLkl,[c_1,c_2],[0.0,0.0],disp=0) x = [d,h, n,reps, i, mle] print x mles.append(x) #mles ︡d1e4d714-2219-438a-98dd-44fcfee14691︡{"stdout":"[2, 6, 127, 1, 0, (-0.3259079515975824, -0.05752026912960705)]"}︡{"stdout":"\n[2, 6, 127, 1, 1, (-0.3259079515975824, -0.05752026912960705)]"}︡{"stdout":"\n"}︡{"stdout":"CPU time: 2.50 s, Wall time: 2.86 s\n"}︡{"done":true}︡ ︠21f4aa9b-ef71-4d92-a2ce-7eb4a7bd8369s︠ #[2, 6, 127, 1, 0, (-0.3259079515975824, -0.05752026912960705)] #[2, 7, 255, 1, 0, (-0.3694074748854663, -0.10600727506430978)] #[2, 8, 511, 1, 0, (-0.39276938517987087, -0.1329531458925401)] #[2, 9, 1023, 1, 0, (-0.4051981012968259, -0.14770363723074506)] #[2, 10, 2047, 1, 0, (-0.41245772002324016, -0.15622461056846243)] ︡e143a023-df2e-428d-999c-16e004340e77︡{"done":true}︡ ︠35b34b06-4ca7-4d98-9f58-1180e7c58a9c︠ #ts=[transmissionProcessTC(graphs.BalancedTree(2,3).to_directed(),0) for _ in range(5)] ts=[transmissionProcessTC(graphs.BalancedTree(3,2).to_directed(),0) for _ in range(5)] ascii_art(ts[0]) # labelled transmission tree is always the "left-branching 3-shark" ignoring labels ︡13cdfa82-3b9e-49d9-b738-7d24bd566583︡{"stdout":" _____________1______________\n / / \n __________2___________ __3___\n / / / / \n _____5______ __4____ _6___ i4\n / / / / / / \ni0 __7___ __8___ i10 12_ i6 \n / / / / / / \n _9___ i7 10_ i11 i1 i5 \n / / / / \n 11_ i9 i3 i12 \n / / \n i2 i8 \n"}︡︡{"done":true} ︠da36f237-d52e-4ee8-862f-365da50cfd44︠ T = ts[4] #view(T) #latex(T) ︡83666bb6-5df7-4683-9953-700642df05bc︡{"done":true}︡ ︠c605c077-f2c2-4cdb-96f1-3bb7176cb20fi︠ %md ## Toroidal Grid Network ︡671ce6f1-02a1-48b0-877d-9a8c43c8cd90︡{"done":true,"md":"## Toroidal Grid Network"} ︠9e0be8ea-799e-4b7f-9a42-f83685fd9058︠ def toroidal2DGrid(n): G = graphs.GridGraph([n,n]) #G.show() # long time G.add_edges([((i,0),(i,n-1)) for i in range(n)]) G.add_edges([((0,i),(n-1,i)) for i in range(n)]) G.relabel() return G ︡f5400699-bd99-451e-8670-b16f4a1b7d57︡{"done":true}︡ ︠105e838e-46f5-4915-afa2-29bad9cd5026︠ ts=[transmissionProcessTC(toroidal2DGrid(10).to_directed(),0) for _ in range(10)] # gives several distinct unlabelled trees depending on the sequence of infection events d=CountsDict([forgetAllLabels(t) for t in ts]) len(d) d.values() ︡514bd13e-ec7c-428a-9748-80e34b2bb03c︡{"stdout":"10\n"}︡{"stdout":"[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n"}︡{"done":true}︡ ︠c25083a2-e3cb-40a6-87c8-8b18a060ecaf︠ d=CountsDict([justTree(t) for t in ts]) # gives several distinct trees depending on the sequence of infection events len(d) d.values() ︡23f99dfe-7b64-4a56-99d8-edde8c9b863d︡{"stdout":"10\n"}︡{"stdout":"[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n"}︡{"done":true}︡ ︠e9c34715-bd10-4387-85e4-a2815873f4c1︠ %time # demo of the mle - way too long... c_1 = lambda p: p[0]+0.9999999 # constraint for alpha > -1 c_2 = lambda p: p[1]+0.9999999 # constraint for beta > -1 # simulation settings m=10 n=m^2 reps=10 repeatMLE=5 for i in range(repeatMLE): ts=[transmissionProcessTC(toroidal2DGrid(m).to_directed(),0) for _ in range(reps)] spc=splitPairsCounts(ts) def negLkl(AB): return negLogLkl_SplitPairCounts2(spc,AB[0],AB[1]) mle=minimize_constrained(negLkl,[c_1,c_2],[0.0,0.0],disp=0) print [i, n,reps,mle] ︡e57916f7-b334-47bd-aa70-7fe553d5445c︡{"stdout":"[0, 100, 10, (-0.7282546466954618, -0.25283759776817866)]"}︡{"stdout":"\n[1, 100, 10, (-0.7292292718997252, -0.2201738282763696)]"}︡{"stdout":"\n[2, 100, 10, (-0.7391465680059263, -0.34600059265553634)]"}︡{"stdout":"\n[3, 100, 10, (-0.7646102471001511, -0.3748667219996737)]"}︡{"stdout":"\n[4, 100, 10, (-0.7249471325940556, -0.20176053588330745)]"}︡{"stdout":"\n"}︡{"stdout":"CPU time: 130.43 s, Wall time: 132.76 s\n"}︡{"done":true}︡ ︠edcbfc88-e4ff-4c7e-a0cc-d6e1c1cd5f34i︠ %md Let us see graph in latex. ︡df2644eb-b83a-4111-9857-4fa0091049ed︡{"done":true,"md":"Let us see graph in latex."} ︠360a99dc-3c52-4dbe-b903-a6d3c12c49f7︠ H = toroidal2DGrid(3).to_directed() H.set_latex_options( graphic_size=(5,5), vertex_size=0.2, edge_thickness=0.04, edge_color='green', vertex_color='green', vertex_label_color='red', format='tkz_graph', #tkz_style='Art', #layout='acyclic' ) ︡f8af87d5-c085-45cc-a2b7-fcfaf7464bd8︡{"done":true}︡ ︠113c223e-6f4f-4f5c-9212-bbbfc1e880e7︠ #view(H) #latex(H) ︡b0789812-417d-4a75-b90e-38e38a045696︡{"done":true}︡ ︠c40efc4a-1151-4f32-939f-6600043ce09b︠ ts=[transmissionProcessTC(toroidal2DGrid(3).to_directed(),0) for _ in range(10)] ︡24b4aa6a-7035-45e8-9323-70d0fc8397ab︡{"done":true}︡ ︠966a9283-f06a-43f4-aa23-9373ac4386d9︠ #view(ts[0]) ︡f7e34156-d9f4-481e-8f03-20e7a407654d︡{"done":true}︡ ︠fff15612-f9e3-4329-b912-ca7fcd40c3d5︠ #latex(ts[0]) ︡efb09ec1-55ab-4195-8e18-f016fd1d25e4︡{"done":true}︡ ︠15aeb500-39fd-4f93-bd8f-bcca326ac8d7︠ #view(ts[1]) ︡020534ad-4863-42c7-b119-216ef6a6785a︡{"done":true}︡ ︠222d6f29-c54e-4cf0-bff9-c2809e98f614︠ #latex(ts[1]) ︡05a8cf71-21ee-4a6b-acda-e4e77d0f97c8︡{"done":true}︡ ︠c92526d8-583e-4290-b432-13e636cab05b︠ #view(ts[3]) #latex(T) ︡729e1b71-0388-436d-bea7-a00b4d1cbf68︡{"done":true}︡ ︠9dd5cc70-ac08-42d9-b453-e6f85d0ef662︠ #latex(ts[3]) ︡2c8d0215-fd37-4224-9906-bf257ffea4a2︡{"done":true}︡ ︠75866b33-324b-471e-ad12-8119ba7f8e11i︠ %md 3D Toroidal Grid ︡b0484e5c-e50c-40f8-8f1b-298464edf435︡{"done":true,"md":"3D Toroidal Grid"} ︠0db73761-b94e-42ad-8835-72f776cc4249︠ def toroidal3DGrid(n): G = graphs.GridGraph([n,n,n]) #G.show() # long time G.add_edges([((0,i,j),(n-1,i,j)) for i in range(n) for j in range(n)]) G.add_edges([((i,0,j),(i,n-1,j)) for i in range(n) for j in range(n)]) G.add_edges([((i,j,0),(i,j,n-1)) for i in range(n) for j in range(n)]) G.relabel() return G ︡9d73e62b-3dae-4461-9de2-0c1025d051ad︡{"done":true}︡ ︠5ceb87b6-df44-47bd-b0f4-e4a1d573127c︠ %time # demo of the mle c_1 = lambda p: p[0]+0.9999999 # constraint for alpha > -1 c_2 = lambda p: p[1]+0.9999999 # constraint for beta > -1 m=10 n=m^3 reps=1 # all reps have the same unlabelled transmission tree topology for i in range(5): ts=[transmissionProcessTC(toroidal3DGrid(m).to_directed(),0) for _ in range(reps)] spc=splitPairsCounts(ts) def negLkl(AB): return negLogLkl_SplitPairCounts2(spc,AB[0],AB[1]) mle=minimize_constrained(negLkl,[c_1,c_2],[0.0,0.0],disp=0) print [n,reps, i, mle] ︡8b6a6586-561f-4887-9de0-a0295bff7197︡{"stdout":"[1000, 1, 0, (-0.6771966619586535, -0.3721134969526218)]"}︡{"stdout":"\n[1000, 1, 1, (-0.6972972063571803, -0.38528953496008167)]"}︡{"stdout":"\n[1000, 1, 2, (-0.7112738575413484, -0.34137924850162177)]"}︡{"stdout":"\n[1000, 1, 3, (-0.6613671977776394, -0.33657729645960566)]"}︡{"stdout":"\n[1000, 1, 4, (-0.695638773816162, -0.34113181129019937)]"}︡{"stdout":"\n"}︡{"stdout":"CPU time: 317.66 s, Wall time: 337.35 s\n"}︡{"done":true}︡ ︠c5f362bd-e330-4f4b-af6a-09e6e1a4d055i︠ %md # Some Random Contact Networks ︡767bc15e-6d4e-4528-9395-88edc5db6248︡{"done":true,"md":"# Some Random Contact Networks"} ︠f33003ca-77c4-4a17-b828-7d73102e0993i︠ %md ## Erdos-Renyi ︡e6d3f622-0686-4e58-9044-961c801cf037︡{"done":true,"md":"## Erdos-Renyi"} ︠9b7844fc-0f09-414a-93a3-d26f7ee8d18es︠ def ErdosReyniConnectedCompOf0(n,p,reps): '''return reps many transmission trees from the connected component containing initial infection vertex 0''' ts=[] i=0; MAXTrials=10000; successfulTrials=0; while (successfulTrials