Indices of events and bSFS dimensions

[A-J-F-Mackintosh]

Hi Gertjan,

I have a few agemo questions and thought it would be best to ask them here (in case they are useful to anyone else).

I followed the example in the docs to set up the IM model and generate an expected bSFS. In setting up the model I defined indices for both the migration event and discrete population-split event.

# populations AB, A, B will have indices 0, 1, 2
sample_configuration = [(), ('a', 'a'), ('b', 'b')]

# index for the migration event
mig_idx = len(sample_configuration)

# index for the split event
split_idx = mig_idx + 1

# put the two events in a list
events = [agemo.MigrationEvent(mig_idx, 1, 2), agemo.PopulationSplitEvent(split_idx, 0, 1, 2)]

My understanding is that these indices need to match how we set up the bSFS evaluator later.

theta = 0.8
theta_along_branchtypes = np.full(len(btc), theta, dtype=np.float64)
params = np.array([1.0, 1.0, 1.0, 0.0])
var = np.hstack([params, theta_along_branchtypes])
bsfs_array = bsfs_eval.evaluate(theta, var, time=2.5)

The above seems to work as expected.

However, I initially thought that I should write params = np.array([C_AB, C_A, C_B, M, T])

Why don't we include T in params given that we have specified an index for it above?

T is instead specified using time. Perhaps bsfs_eval already knows the index of the discrete event?

If I had set the split index as 3 and the migration index as 4, then would it be possible to set up the model, i.e. how would I specify the params?

I tried a few different parameter combinations and I think that the returned bSFS has dimensions [hetA, hetB, fixed, hetAB], which differs from the representation that gIMble uses (at least back in my day) - [hetB, hetA, hetAB, fixed]. Is that correct?

Many thanks and best wishes,

Alex

[GertjanBisschop]

Hi @A-J-F-Mackintosh,
Yep, as always a couple of very good questions!
Firstly, the easiest question to answer. Yes, the default order of branch types is different from the order used in gimble. This is because in agemo all branch types are represented as a unique sequence of 0's and 1's. Their final order in the bSFS vector is determined by a simple sorting of those representations. This is not how we came up with the order in gimble.
However, you can specify your own order by providing a branchtype_dict that maps each branchtype to a unique integer index when making a BranchTypeCounter().

Now onto the issue of indices and parameters. I agree that this is currently a bit awkward. It should be possible though to set the split index to 3, and the migration index to 4. In the future, indices should be set automatically based on the order in which you insert them. I am still a bit on the fence as to what the best way would be to pass the parameter values once you have specified a model.

Let me know if there are other things that are slightly annoying. I am hoping to polish off some of the rougher edges in the near future.

[A-J-F-Mackintosh]

Hi Gertjan,

Really helpful answers, thank you. I agree that setting the indices automatically would be a good way to avoid any confusion.

A follow-up question - I am interested in setting up a model where lineages are sampled from a single population and then there is a discrete event at time T where N_e changes. I think this is possible to set up by adjusting the IM model (setting M to zero and only sampling lineages from one population) but this seems like quite an awkward way of doing it.

Does agemo have a class for that kind of event yet? I had a look in the docs and code but couldn't see what I was looking for.

I realise that this is tangential to the original issue, so happy to open a new one or chat by email.

Thanks again,

Alex

[GertjanBisschop]

I think this can be achieved by means of the PopulationSplitEvent class where derived only contains a single population.

sample_configuration = [(, ), ('a', 'a')]
btc = agemo.BranchTypeCounter(sample_configuration)
idx = len(sample_configuration)
events_list = [agemo.PopulationSplitEvent(idx, 0, 1)]

It would probably make more sense to implement a general MassMigrationEvent class instead that allows you to move one or more populations into another one. This can then be used by PopulationSplitEvent. Feel free to have a stab at this. This should be straightforward. If you're up for this, don't hesitate to ask for pointers. I'll open an issue to track this.

[A-J-F-Mackintosh]

Using PopulationSplitEvent like that seems to work, thanks!

I think adding a MassMigrationEvent would be useful and am happy to give it a go.

I think the main pointer I would need is about the method _single_step() that all events have. My understanding is that this must return a tuple of (idx, count, state) where idx is the index of the event and state is the sample configuration after the event. I am not sure about count though and why it is needed for the multiplier_array.

I would also be interested in making a class for an exponentially distributed event that results in a burst of coalescence (emulating a recurrent sweep process). Perhaps this could already be done with CoalescenceEvent() though? I did try, but ran into problems because the _single_step() method for this class does not return an index. Should be possible though, right?

[GertjanBisschop]

Right, so you would want a MultiMergerCoalescenceEvent essentially, right? Should be straightforward indeed. We can add an argument that specifies the number of lineages that can coalesce. Regular coalescence events would have a parameter equalling 2. I'll open another issue on that one. It would be good to coordinate with @DerekSetter. I think we'll need that functionality for the sweeps anyway.
Coalescence events are further grouped in the CoalescenceEventsSuite class. Did you look into that?