evaluate msprime's input reqs; discuss algos to produce CRs #1

ashander · 2016-11-10T00:51:36Z

No description provided.

ashander

Discussion points for today

ashander · 2016-11-10T00:53:48Z

devel/test_msprime.py

+
+### Works!!  Compare pictures of trees.
+# Node identification of internal nodes not the same though:
+


ashander · 2016-11-10T00:55:17Z

devel/test_msprime.py

+    msprime.CoalescenceRecord(left=0.2,  right=0.8,  node=4,  children=(0,  3),  time=0.2,  population=0),
+    msprime.CoalescenceRecord(left=0.0,  right=0.2,  node=5,  children=(2,  4),  time=0.4,  population=0),
+    msprime.CoalescenceRecord(left=0.2,  right=0.8,  node=5,  children=(0,  2),  time=0.4,  population=0),
+    msprime.CoalescenceRecord(left=0.8,  right=1.0,  node=5,  children=(2,  4),  time=0.4,  population=0),


I don't fully understand what's happening to node 4 here

ashander · 2016-11-10T17:00:33Z

devel/test_msprime.py

+    msprime.CoalescenceRecord(left=0.8,  right=1.0,  node=5,  children=(1,  2),  time=0.5,  population=0),
+    msprime.CoalescenceRecord(left=0.8,  right=1.0,  node=6,  children=(0,  5),  time=0.7,  population=0),
+    msprime.CoalescenceRecord(left=0.0,  right=0.2,  node=7,  children=(0,  5),  time=1.0,  population=0),
+    ]


I used this example in reading through Generating trees. The ordered index vectors would be (indexing starting at 1)
I = [1, 4, 8, 2, 5, 3, 6, 7] and O = [8, 4, 1, 5, 2, 7, 6, 3].

initialize by setting pi_j = 0 for j in (0, 7). So `pi = [0, 0 ,0, 0, 0, 0, 0, 0]

Insert first three records for h = 1, 4, 8

h=1, pi = [0, 0, 4, 4, ...]

h=4, pi = [0, 5, 5, 4, ...] # note we reset 2 here!

h =8, pi = [7, 5, 5, 4, 0, 7, ..]

visit the tree -- but 4 isn't properly encoded!

Replace all node 5 with: CoalescenceRecords(0, 1.0, 5, (1,4), 0.5) at index j=4.

problem was no records overlapping 0.0 has a parent for 4

should be: (?)

my_records = [ msprime.CoalescenceRecord(left=0.0, right=0.2, node=4, children=(2, 3), time=0.4, population=0), # 1 msprime.CoalescenceRecord(left=0.2, right=0.8, node=4, children=(0, 2), time=0.4, population=0), # 2 msprime.CoalescenceRecord(left=0.8, right=1.0, node=4, children=(2, 3), time=0.4, population=0), # 3 msprime.CoalescenceRecord(left=0.0, right=1.0, node=5, children=(1, 4), time=0.5, population=0), # 4 msprime.CoalescenceRecord(left=0.8, right=1.0, node=6, children=(0, 5), time=0.7, population=0), # 5 msprime.CoalescenceRecord(left=0.0, right=0.2, node=7, children=(0, 5), time=1.0, population=0) # 6 ]

Then indices 7 and 8 become 4 and 5 so I = [1, 4, 6, 2, 3, 5] and O = [ 6, 1, 2, 5, 4, 3]

ashander · 2016-11-10T17:28:28Z

devel/test_msprime.py

+    msprime.CoalescenceRecord(left=0.8,  right=1.0,  node=4,  children=(1,  2),  time=0.5,  population=0),
+    msprime.CoalescenceRecord(left=0.8,  right=1.0,  node=5,  children=(0,  4),  time=0.7,  population=0),
+    msprime.CoalescenceRecord(left=0.0,  right=0.2,  node=6,  children=(0,  4),  time=1.0,  population=0)]
+


Lets apply alg T:
The ordered index vectors would be (indexing starting at 1)
I = [2, 6, 1, 3, 4, 5] and O = [6, 2, 3, 1, 5, 4].

initialize by setting pi_j = 0 for j in 1 to 6 + 1. So pi = [0, 0 ,0, 0, 0, 0, 0] and j = 1, k = 1, x= 0.0. Also M = 6

Insert for x =0.0?

j = 1; h=2, insert pi = [0, 4 ,4, 0, 0, 0, 0]

j = 2; h=6, insert pi = [6, 4 ,4, 0, 6, 0, 0]

j = 3; h=1; but x != l_1 goto visit

visit the tree [6, 4 ,4, 0, 6, 0, 0] -- correct but no root marked! set x = l_1 = 0.2

remove for x=0.2?

k = 1; h = 6; because r_6 == x we remove: [0, 4 ,4, 0, 0, 0, 0]

k = 2; h =2; because r_2 == x we remove: [0, 0, 0, 0, 0, 0, 0]

k = 3; h = 3; because r_3 != x goto Insert.

insert for x=0.2?

j = 3; h = 1. pi = [3, 0, 3, 0, 0, 0, 0]

j = 4; h = 3. pi = [3, 4, 3, 4, 0, 0, 0]

j = 5; h = 4 but x != l_4 = 0.8. goto visit

visit the tree [3, 4, 3, 4, 0, 0, 0]-- correct but no root marked! set x = l_4 = 0.8

remove for x=0.8?

k = 3; h = 3; because r_3 == x we remove: [3, 0, 3, 0, 0, 0, 0]

k = 4; h = 1; because r_1 == x we remove: [0, 0, 0, 0, 0, 0, 0]

k = 5; h = 5; because r_5 != x goto Insert

insert for x=0.8?

j = 5; h = 4. pi = [0, 4, 4, 0, 0, 0, 0]

j = 6; h = 5. pi = [5, 4, 4, 5, 0, 0, 0]

j = 7

visit [5, 4, 4, 5, 0, 0, 0] --- correct (but no root). because j > M terminate.

Only a single rec for (1,4) -> 5 on (0.0, 1.0)

ashander · 2016-11-10T19:37:01Z

OK I've applied T to the corrected records an it's unclear why it would fail. This is currently in devel/notes.md. I'll try it in code and report back

The output from trees py is in notes.md and lists the records in and out during alg T. It matches with what I did by hand The reasons are described in notes.md

ashander · 2016-11-10T21:04:44Z

devel/notes.md

+children (beyond the first child) added or removed is
+equal.
+For binary trees (like we have) this is just that the number of records in and
+out match.


see line 2225 (my comment) of this diff in msprime for maybe some more detail on the underlying reasons ashander/msprime@1fcfb2e081cb34ac1ef4092a056b4c407f34c1fbx

link doesn't work. very nice writeup tho

So: For any rooted tree, the sum over nodes of (the number of children minus one)
is equal to the number of tips minus one. (proof: order the tips, then from each draw the path from that tip towards the root until it hits an already-drawn path; associate each tip with the last edge in its path.)
If we've just constructed the first tree, in_count is this number.

This suggests that the above tree sequence should fail the first_tree check, since self-> sample_size should be 3, while the tree has 4 tips.

I've run with the more informative errors, and contrary to this comment yesterday, the first_tree check passes for these records. Failure is still at 12. :?

ashander · 2016-11-10T22:23:21Z

oh the one to the diff in msprime? oops. I'm at the vet.
On Thursday, November 10, 2016, Peter Ralph notifications@github.com
wrote:

@petrelharp commented on this pull request.

In devel/notes.md #1:

} else {

if (in_count != out_count) {

ret = MSP_ERR_BAD_COALESCENCE_RECORDS_12;

goto out;

}

+The in_count and
+out_count
+both accumulate the number of children (less 1 for some reason) in or out when adding or removing a set of records.
+
+The above snippet appears just before 'returning' the tree and seems to be
+a check that, at any given point after the first tree, the cumulative number of
+children (beyond the first child) added or removed is
+equal.
+For binary trees (like we have) this is just that the number of records in and
+out match.

link doesn't work. very nice writeup tho

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
#1, or mute the thread
https://github.com/notifications/unsubscribe-auth/AAfLOAXZCYcRbX3x7AHvyZCl2r9PjVMuks5q85W9gaJpZM4KuKBq
.

-Jaime

ashander · 2016-11-10T22:35:47Z

I can't fully respond right now but maybe this link works. look for my
comment in tree_sequence.c
ashander/msprime@dc8b8cc

On Thursday, November 10, 2016, Peter Ralph notifications@github.com
wrote:

@petrelharp commented on this pull request.

In devel/notes.md #1:

} else {

if (in_count != out_count) {

ret = MSP_ERR_BAD_COALESCENCE_RECORDS_12;

goto out;

}

+The in_count and
+out_count
+both accumulate the number of children (less 1 for some reason) in or out when adding or removing a set of records.
+
+The above snippet appears just before 'returning' the tree and seems to be
+a check that, at any given point after the first tree, the cumulative number of
+children (beyond the first child) added or removed is
+equal.
+For binary trees (like we have) this is just that the number of records in and
+out match.

So: For any rooted tree, the sum over nodes of (the number of children
minus one)
is equal to the number of tips minus one. (proof: order the tips, then
from each draw the path from that tip towards the root until it hits an
already-drawn path; associate each tip with the last edge in its path.)
If we've just constructed the first tree, in_count is this number.

This suggests that the above tree sequence should fail the first_tree
check, since self-> sample_size should be 3, while the tree has 4 tips.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
#1, or mute the thread
https://github.com/notifications/unsubscribe-auth/AAfLOCd28SZfb5HbGHsEO7KcncmtcEtaks5q85rDgaJpZM4KuKBq
.

-Jaime

also, add results with error located

msprime infers 3 is a root in these records

in 'adding whole genome info'

ashander · 2016-11-11T19:58:25Z

devel/test_msprime.py

+#  	in: CoalescenceRecord(left=0.8, right=1.0, node=5, children=(2, 4), time=0.4, population=0)
+#  	in: CoalescenceRecord(left=0.8, right=1.0, node=7, children=(0, 6), time=0.7, population=0)
+#  ([7, 6, 5, -1, 5, 6, 7, -1, -1], [[], [], [], [], [], (2, 4), (1, 5), (0, 6), []])
+


Incorrect trees inferred here -- 3 is not child to 5 -- need to check the CRs are as intended

ashander · 2016-11-12T00:54:50Z

devel/algorithm-notes.md

+
+```
+
+Note that we've output the records not quite in time order.


if I understand correctly ( @petrelharp let me know if not) part of the reason for this, and for not assigning new IDs immediately after a split as in Jerome's outline (which I pasted in below) is that we want lineages to persist across multiple generations

or is it for the purpose of good renumbering at the end? unsure as of now

oh, good point. that would work. I was trying to be a bit too clever, and getting tripped up by the idea of individuals existing across multiple linages.

So: we could assign new IDs right after a split. It seems like they'd be able to persist across multiple gens fine (eg as in his example for individual ID 1). that might come at the cost of needing to be more careful at the end when we renumber but unclear it will be a big deal

ashander · 2016-11-12T01:22:25Z

devel/algorithm-notes.md

+
+... where we only really have to output the coalescence record 
+so as to remember that lineage 9 is a continuation of lineage 3.
+(But, either works.)


Either works, but in this instance because b died so 8 is going to be 'dangling' tip that isn't a sample

which is fine

Yes. OK one Q: if we weren't tracking individuals (in tree) but rather just lineages by ID number what would we lose?

not sure what you mean

interactive noodling

3547fc7

ashander commented Nov 10, 2016

View reviewed changes

ashander and others added 3 commits November 10, 2016 10:57

Only a single rec for (1,4) -> 5 on (0.0, 1.0)

7a7f414

Alg T applied. Unclear why it'd fail :?

0243c66

Merge pull request #2 from ashander/fix/plr-records

838a72c

Only a single rec for (1,4) -> 5 on (0.0, 1.0)

ashander added 2 commits November 10, 2016 11:41

Jerome's trees.py

0f37dc1

Works in python; fails in c; I think we know why :D

98e34da

The output from trees py is in notes.md and lists the records in and out during alg T. It matches with what I did by hand The reasons are described in notes.md

ashander commented Nov 10, 2016

View reviewed changes

petrelharp added 2 commits November 11, 2016 10:11

python algorithm works with partial dangling tips;

a204580

notes on algorithms

53001e1

ashander changed the title ~~discuss interactive noodling~~ evaluate msprime's input reqs; discuss algos to produce CRs Nov 11, 2016

ashander added 3 commits November 11, 2016 11:36

Final set of records still work in C

40ecd70

also, add results with error located

Why 'adding whole genome info for phantom' fails

7d1e668

msprime infers 3 is a root in these records

Alg gives incorrect trees for the records

c2fa7fa

in 'adding whole genome info'

ashander commented Nov 11, 2016

View reviewed changes

petrelharp and others added 4 commits November 11, 2016 15:18

notes oin time

2aa174c

Close some math envs

1bee243

example

b68456c

Add Jerome's email outlining fwd time

4100b16

ashander commented Nov 12, 2016

View reviewed changes

ashander and others added 3 commits November 11, 2016 17:15

Close a math env

570d3cf

modified example

54bc106

modification to bookkeeping

7ac5c33

ashander commented Nov 12, 2016

View reviewed changes

petrelharp added 2 commits November 11, 2016 17:31

algorithm, no recomb

4f233e9

python tree iterator works with single-offspring coalescence records

239ddac

petrelharp added 2 commits November 11, 2016 17:59

.

c568ec6

an algorithm, that probably works

1e0360f

petrelharp merged commit 1e0360f into master Nov 13, 2016

petrelharp deleted the plr branch November 18, 2016 17:48


		### Works!! Compare pictures of trees.
		# Node identification of internal nodes not the same though:


		```

		Note that we've output the records not quite in time order.

evaluate msprime's input reqs; discuss algos to produce CRs #1

evaluate msprime's input reqs; discuss algos to produce CRs #1

Conversation

ashander commented Nov 10, 2016

ashander left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ashander Nov 10, 2016 • edited Loading

Choose a reason for hiding this comment

ashander Nov 10, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

petrelharp Nov 10, 2016 • edited by ashander Loading

Choose a reason for hiding this comment

ashander Nov 10, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ashander commented Nov 10, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ashander commented Nov 10, 2016

@petrelharp commented on this pull request.

} else {

if (in_count != out_count) {

ret = MSP_ERR_BAD_COALESCENCE_RECORDS_12;

goto out;

}

ashander commented Nov 10, 2016

@petrelharp commented on this pull request.

} else {

if (in_count != out_count) {

ret = MSP_ERR_BAD_COALESCENCE_RECORDS_12;

goto out;

}

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ashander Nov 10, 2016 •

edited

Loading

ashander Nov 10, 2016 •

edited

Loading

petrelharp Nov 10, 2016 •

edited by ashander

Loading

ashander Nov 10, 2016 •

edited

Loading