Indices of events and bSFS dimensions #18

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